PWGCF/Correlations/Base/AliUEHist.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliUEHist.cxx 20164 2007-08-14 15:31:50Z morsch $ */
  17
  18 //
  19 //
  20 // encapsulate histogram and corrections for one underlying event histogram
  21 //
  22 //
  23 // Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero
  24
  25 #include "AliUEHist.h"
  26 #include "AliCFContainer.h"
  27 #include "THnSparse.h"
  28 #include "TMath.h"
  29 #include "TList.h"
  30 #include "TCollection.h"
  31 #include "TH1D.h"
  32 #include "TH2D.h"
  33 #include "TH3D.h"
  34 #include "AliLog.h"
  35 #include "TCanvas.h"
  36 #include "TF1.h"
  37 #include "AliTHn.h"
  38 #include "THn.h"
  39
  40 ClassImp(AliUEHist)
  41
  42 const Int_t AliUEHist::fgkCFSteps = 11;
  43
  44 AliUEHist::AliUEHist(const char* reqHist) :
  45   TObject(),
  46   fkRegions(4),
  47   fEventHist(0),
  48   fTrackHistEfficiency(0),
  49   fFakePt(0),
  50   fEtaMin(0),
  51   fEtaMax(0),
  52   fPtMin(0),
  53   fPtMax(0),
  54   fCentralityMin(0),
  55   fCentralityMax(0),
  56   fZVtxMin(0),
  57   fZVtxMax(0),
  58   fContaminationEnhancement(0),
  59   fCombineMinMax(0),
  60   fTrackEtaCut(0),
  61   fCache(0),
  62   fGetMultCacheOn(kFALSE),
  63   fGetMultCache(0),
  64   fHistogramType(reqHist)
  65 {
  66   // Constructor
  67
  68   for (UInt_t i=0; i<fkRegions; i++)
  69     fTrackHist[i] = 0;
  70
  71   if (strlen(reqHist) == 0)
  72     return;
  73
  74   Printf("Creating AliUEHist with %s", reqHist);
  75
  76   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  77   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  78
  79   const char* title = "";
  80
  81   // track level
  82   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  83   Int_t iTrackBin[6];
  84   Double_t* trackBins[6];
  85   const char* trackAxisTitle[6];
  86
  87   // eta
  88   const Int_t kNEtaBins = 20;
  89   Double_t etaBins[20+1];
  90   for (Int_t i=0; i<=kNEtaBins; i++)
  91     etaBins[i] = -1.0 + 0.1 * i;
  92   const char* etaTitle = "#eta";
  93   iTrackBin[0] = kNEtaBins;
  94   trackBins[0] = etaBins;
  95   trackAxisTitle[0] = etaTitle;
  96
  97   // delta eta
  98   const Int_t kNDeltaEtaBins = 48;
  99   Double_t deltaEtaBins[kNDeltaEtaBins+1] = { -2.4, -2.3, -2.2, -2.1,
 100                                               -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1,
 101                                               -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
 102                                               0,
 103                                               0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
 104                                               1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
 105                                               2.1, 2.2, 2.3, 2.4 };
 106
 107   const Int_t kNDeltaEtaBinsTTR = 40+4;
 108   Double_t deltaEtaBinsTTR[kNDeltaEtaBinsTTR+1] = { -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1,
 109                                                   -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
 110                                                   -0.05, -0.025, 0, 0.025, 0.05,
 111                                                   0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
 112                                                   1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0 };
 113
 114   // pT
 115   iTrackBin[1] = 9;
 116   Double_t pTBins[] = {0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0 };
 117   trackBins[1] = pTBins;
 118   trackAxisTitle[1] = "p_{T} (GeV/c)";
 119
 120   // pT, fine
 121   const Int_t kNpTBinsFine = 22;
 122   Double_t pTBinsFine[kNpTBinsFine+1] = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0};
 123
 124   // pT, wide range
 125   //const Int_t kNpTBinsWideRange = 20;
 126   //Double_t pTBinsWideRange[kNpTBinsWideRange+1] = {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, 10.0, 12.0, 15.0};
 127
 128   // pT,lead binning 1
 129   const Int_t kNLeadingpTBins = 100;
 130   Double_t leadingpTBins[kNLeadingpTBins+1];
 131   for (Int_t i=0; i<=kNLeadingpTBins; i++)
 132     leadingpTBins[i] = 0.5 * i;
 133
 134   // pT,lead binning 2
 135   //const Int_t kNLeadingpTBins2 = 8;
 136 //   Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0 };
 137   const Int_t kNLeadingpTBins2 = 6;
 138   Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0 };
 139
 140   // phi,lead; this binning starts at -pi/2 and is modulo 3
 141   const Int_t kNLeadingPhiSpacing = 72;
 142   const Int_t kNLeadingPhiBins = kNLeadingPhiSpacing;
 143   Double_t leadingPhiBins[kNLeadingPhiBins+1];
 144   for (Int_t i=0; i<=kNLeadingPhiSpacing; i++)
 145     leadingPhiBins[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * i;
 146
 147   const Int_t kNLeadingPhiBinsTTR = kNLeadingPhiSpacing+4;
 148   Double_t leadingPhiBinsTTR[kNLeadingPhiBinsTTR+1];
 149   for (Int_t i=0; i<=17; i++)
 150     leadingPhiBinsTTR[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * i;
 151   leadingPhiBinsTTR[18] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2;
 152   leadingPhiBinsTTR[19] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 1.5;
 153   leadingPhiBinsTTR[20] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 2; // = 0
 154   leadingPhiBinsTTR[21] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 2.5;
 155   leadingPhiBinsTTR[22] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 3;
 156   for (Int_t i=19; i<=kNLeadingPhiSpacing; i++)
 157     leadingPhiBinsTTR[i+4] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * i;
 158
 159   // multiplicity
 160   const Int_t kNMultiplicityBins = 15;
 161   Double_t multiplicityBins[kNMultiplicityBins+1];
 162   for (Int_t i=0; i<=kNMultiplicityBins; i++)
 163     multiplicityBins[i] = -0.5 + i;
 164   multiplicityBins[kNMultiplicityBins] = 200;
 165
 166   trackBins[3] = multiplicityBins;
 167   iTrackBin[3] = kNMultiplicityBins;
 168   trackAxisTitle[3] = "multiplicity";
 169
 170   // centrality (in %)
 171   const Int_t kNCentralityBins = 5 + 1 + 9;
 172   Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
 173
 174   const Int_t kNCentralityBinsCourse = 5;
 175   Double_t centralityBinsCourse[] = { 0, 20, 40, 60, 80, 100.1 };
 176
 177   // particle species
 178   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
 179   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
 180
 181   // vtx-z axis
 182   const char* vertexTitle = "z-vtx (cm)";
 183   const Int_t kNVertexBins = 7;
 184   Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
 185   const Int_t kNVertexBins2 = 10;
 186   Double_t vertexBins2[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10 };
 187
 188   Int_t useVtxAxis = 0;
 189   Bool_t useTTRBinning = kFALSE;
 190   Bool_t useCourseCentralityBinning = kFALSE;
 191
 192   // selection depending on requested histogram
 193   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
 194   if (strcmp(reqHist, "NumberDensitypT") == 0)
 195   {
 196     axis = 0;
 197     title = "d^{2}N_{ch}/d#varphid#eta";
 198   }
 199   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 200   {
 201     axis = 1;
 202     title = "d^{2}N_{ch}/d#varphid#eta";
 203   }
 204   else if (TString(reqHist).BeginsWith("NumberDensityPhiCentrality"))
 205   {
 206     if (TString(reqHist).Contains("Vtx"))
 207       useVtxAxis = 1;
 208     if (TString(reqHist).Contains("Vtx10"))
 209       useVtxAxis = 2;
 210     if (TString(reqHist).Contains("TTR"))
 211       useTTRBinning = kTRUE;
 212     if (TString(reqHist).Contains("Course"))
 213       useCourseCentralityBinning = kTRUE;
 214
 215     reqHist = "NumberDensityPhiCentrality";
 216     fHistogramType = reqHist;
 217     axis = 2;
 218     title = "d^{2}N_{ch}/d#varphid#eta";
 219   }
 220   else if (strcmp(reqHist, "SumpT") == 0)
 221   {
 222     axis = 0;
 223     title = "d^{2}#Sigma p_{T}/d#varphid#eta";
 224   }
 225   else
 226     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 227
 228   UInt_t initRegions = fkRegions;
 229
 230   if (axis == 0)
 231   {
 232     trackBins[2] = leadingpTBins;
 233     iTrackBin[2] = kNLeadingpTBins;
 234     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 235
 236   }
 237   else if (axis == 1)
 238   {
 239     nTrackVars = 5;
 240     initRegions = 1;
 241
 242     iTrackBin[2] = kNLeadingpTBins2;
 243     trackBins[2] = leadingpTBins2;
 244     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 245
 246     iTrackBin[4] = kNLeadingPhiBins;
 247     trackBins[4] = leadingPhiBins;
 248     trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
 249   }
 250   else if (axis == 2)
 251   {
 252     nTrackVars = 5;
 253     initRegions = 1;
 254
 255     iTrackBin[0] = (useTTRBinning) ? kNDeltaEtaBinsTTR : kNDeltaEtaBins;
 256     trackBins[0] = (useTTRBinning) ? deltaEtaBinsTTR : deltaEtaBins;
 257     trackAxisTitle[0] = "#Delta#eta";
 258
 259     iTrackBin[2] = kNLeadingpTBins2;
 260     trackBins[2] = leadingpTBins2;
 261     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 262
 263     trackBins[3] = (useCourseCentralityBinning) ? centralityBinsCourse : centralityBins;
 264     iTrackBin[3] = (useCourseCentralityBinning) ? kNCentralityBinsCourse : kNCentralityBins;
 265     trackAxisTitle[3] = "centrality";
 266
 267     iTrackBin[4] = (useTTRBinning) ? kNLeadingPhiBinsTTR : kNLeadingPhiBins;
 268     trackBins[4] = (useTTRBinning) ? leadingPhiBinsTTR : leadingPhiBins;
 269     trackAxisTitle[4] = "#Delta#varphi (rad.)";
 270
 271     if (useVtxAxis > 0)
 272     {
 273       nTrackVars = 6;
 274       iTrackBin[5] = (useVtxAxis == 1) ? kNVertexBins : kNVertexBins2;
 275       trackBins[5] = (useVtxAxis == 1) ? vertexBins : vertexBins2;
 276       trackAxisTitle[5] = vertexTitle;
 277     }
 278   }
 279
 280   for (UInt_t i=0; i<initRegions; i++)
 281   {
 282     if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
 283       fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 284     else
 285       fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 286
 287     for (Int_t j=0; j<nTrackVars; j++)
 288     {
 289       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 290       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 291     }
 292
 293     SetStepNames(fTrackHist[i]);
 294   }
 295
 296   // event level
 297   Int_t nEventVars = 2;
 298   Int_t iEventBin[3];
 299
 300   // track 3rd and 4th axis --> event 1st and 2nd axis
 301   iEventBin[0] = iTrackBin[2];
 302   iEventBin[1] = iTrackBin[3];
 303
 304   // plus track 5th axis (in certain cases)
 305   if (axis == 2 && useVtxAxis)
 306   {
 307     nEventVars = 3;
 308     iEventBin[2] = iTrackBin[5];
 309   }
 310
 311   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
 312
 313   fEventHist->SetBinLimits(0, trackBins[2]);
 314   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 315
 316   fEventHist->SetBinLimits(1, trackBins[3]);
 317   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 318
 319   if (axis == 2 && useVtxAxis)
 320   {
 321     fEventHist->SetBinLimits(2, trackBins[5]);
 322     fEventHist->SetVarTitle(2, trackAxisTitle[5]);
 323   }
 324
 325   SetStepNames(fEventHist);
 326
 327   iTrackBin[1] = kNpTBinsFine;
 328   iTrackBin[2] = kNSpeciesBins;
 329   iTrackBin[4] = kNVertexBins2;
 330
 331   fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 4, 5, iTrackBin);
 332   fTrackHistEfficiency->SetBinLimits(0, trackBins[0]);
 333   fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]);
 334   fTrackHistEfficiency->SetBinLimits(1, pTBinsFine);
 335   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
 336   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
 337   fTrackHistEfficiency->SetVarTitle(2, "particle species");
 338   fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
 339   fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
 340   fTrackHistEfficiency->SetBinLimits(4, vertexBins2);
 341   fTrackHistEfficiency->SetVarTitle(4, vertexTitle);
 342
 343   fFakePt = new TH3F("fFakePt","fFakePt;p_{T,rec};p_{T};centrality", 200, 0, 20, 200, 0, 20, 20, 0, 100);
 344 }
 345
 346 //_____________________________________________________________________________
 347 AliUEHist::AliUEHist(const AliUEHist &c) :
 348   TObject(),
 349   fkRegions(4),
 350   fEventHist(0),
 351   fTrackHistEfficiency(0),
 352   fFakePt(0),
 353   fEtaMin(0),
 354   fEtaMax(0),
 355   fPtMin(0),
 356   fPtMax(0),
 357   fCentralityMin(0),
 358   fCentralityMax(0),
 359   fZVtxMin(0),
 360   fZVtxMax(0),
 361   fContaminationEnhancement(0),
 362   fCombineMinMax(0),
 363   fTrackEtaCut(0),
 364   fCache(0),
 365   fGetMultCacheOn(kFALSE),
 366   fGetMultCache(0),
 367   fHistogramType()
 368 {
 369   //
 370   // AliUEHist copy constructor
 371   //
 372
 373   ((AliUEHist &) c).Copy(*this);
 374 }
 375
 376 //____________________________________________________________________
 377 void AliUEHist::SetStepNames(AliCFContainer* container)
 378 {
 379   // sets the names of the correction steps
 380
 381   for (Int_t i=0; i<fgkCFSteps; i++)
 382     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 383 }
 384
 385 //____________________________________________________________________
 386 AliUEHist::~AliUEHist()
 387 {
 388   // Destructor
 389
 390   for (UInt_t i=0; i<fkRegions; i++)
 391   {
 392     if (fTrackHist[i])
 393     {
 394       delete fTrackHist[i];
 395       fTrackHist[i] = 0;
 396     }
 397   }
 398
 399   if (fEventHist)
 400   {
 401     delete fEventHist;
 402     fEventHist = 0;
 403   }
 404
 405   if (fTrackHistEfficiency)
 406   {
 407     delete fTrackHistEfficiency;
 408     fTrackHistEfficiency = 0;
 409   }
 410
 411   if (fFakePt)
 412   {
 413     delete fFakePt;
 414     fFakePt = 0;
 415   }
 416
 417   if (fCache)
 418   {
 419     delete fCache;
 420     fCache = 0;
 421   }
 422 }
 423
 424 //____________________________________________________________________
 425 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 426 {
 427   // assigment operator
 428
 429   if (this != &c)
 430     ((AliUEHist &) c).Copy(*this);
 431
 432   return *this;
 433 }
 434
 435 //____________________________________________________________________
 436 void AliUEHist::Copy(TObject& c) const
 437 {
 438   // copy function
 439
 440   AliUEHist& target = (AliUEHist &) c;
 441
 442   for (UInt_t i=0; i<fkRegions; i++)
 443     if (fTrackHist[i])
 444       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 445
 446   if (fEventHist)
 447     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 448
 449   if (fTrackHistEfficiency)
 450     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
 451
 452   if (fFakePt)
 453     target.fFakePt = dynamic_cast<TH3F*> (fFakePt->Clone());
 454
 455   target.fEtaMin = fEtaMin;
 456   target.fEtaMax = fEtaMax;
 457   target.fPtMin = fPtMin;
 458   target.fPtMax = fPtMax;
 459   target.fCentralityMin = fCentralityMin;
 460   target.fCentralityMax = fCentralityMax;
 461   target.fZVtxMin = fZVtxMin;
 462   target.fZVtxMax = fZVtxMax;
 463
 464   if (fContaminationEnhancement)
 465     target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
 466
 467   target.fCombineMinMax = fCombineMinMax;
 468   target.fTrackEtaCut = fTrackEtaCut;
 469   target.fHistogramType = fHistogramType;
 470 }
 471
 472 //____________________________________________________________________
 473 Long64_t AliUEHist::Merge(TCollection* list)
 474 {
 475   // Merge a list of AliUEHist objects with this (needed for
 476   // PROOF).
 477   // Returns the number of merged objects (including this).
 478
 479   if (!list)
 480     return 0;
 481
 482   if (list->IsEmpty())
 483     return 1;
 484
 485   TIterator* iter = list->MakeIterator();
 486   TObject* obj;
 487
 488   // collections of objects
 489   const UInt_t kMaxLists = fkRegions+3;
 490   TList** lists = new TList*[kMaxLists];
 491
 492   for (UInt_t i=0; i<kMaxLists; i++)
 493     lists[i] = new TList;
 494
 495   Int_t count = 0;
 496   while ((obj = iter->Next())) {
 497
 498     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 499     if (entry == 0)
 500       continue;
 501
 502     for (UInt_t i=0; i<fkRegions; i++)
 503       if (entry->fTrackHist[i])
 504         lists[i]->Add(entry->fTrackHist[i]);
 505
 506     lists[fkRegions]->Add(entry->fEventHist);
 507     lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
 508     if (entry->fFakePt)
 509       lists[fkRegions+2]->Add(entry->fFakePt);
 510
 511     count++;
 512   }
 513   for (UInt_t i=0; i<fkRegions; i++)
 514     if (fTrackHist[i])
 515       fTrackHist[i]->Merge(lists[i]);
 516
 517   fEventHist->Merge(lists[fkRegions]);
 518   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
 519   if (fFakePt)
 520     fFakePt->Merge(lists[fkRegions+2]);
 521
 522   for (UInt_t i=0; i<kMaxLists; i++)
 523     delete lists[i];
 524
 525   delete[] lists;
 526
 527   return count+1;
 528 }
 529
 530 //____________________________________________________________________
 531 void AliUEHist::SetBinLimits(THnBase* grid)
 532 {
 533   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 534
 535   if (fEtaMax > fEtaMin)
 536     grid->GetAxis(0)->SetRangeUser(fEtaMin, fEtaMax);
 537   if (fPtMax > fPtMin)
 538     grid->GetAxis(1)->SetRangeUser(fPtMin, fPtMax);
 539 }
 540
 541 //____________________________________________________________________
 542 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 543 {
 544   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 545
 546   SetBinLimits(grid->GetGrid());
 547 }
 548
 549 //____________________________________________________________________
 550 void AliUEHist::ResetBinLimits(THnBase* grid)
 551 {
 552   // resets all bin limits
 553
 554   for (Int_t i=0; i<grid->GetNdimensions(); i++)
 555     if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange))
 556       grid->GetAxis(i)->SetRangeUser(0, -1);
 557 }
 558
 559 //____________________________________________________________________
 560 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 561 {
 562   // resets all bin limits
 563
 564   ResetBinLimits(grid->GetGrid());
 565 }
 566
 567 //____________________________________________________________________
 568 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 569 {
 570   // prints the number of empty bins in the track end event histograms in the given step
 571
 572   Int_t binBegin[4];
 573   Int_t binEnd[4];
 574
 575   for (Int_t i=0; i<4; i++)
 576   {
 577     binBegin[i] = 1;
 578     binEnd[i]   = fTrackHist[0]->GetNBins(i);
 579   }
 580
 581   if (fEtaMax > fEtaMin)
 582   {
 583     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
 584     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
 585   }
 586
 587   if (fPtMax > fPtMin)
 588   {
 589     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
 590     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
 591   }
 592
 593   if (ptLeadMax > ptLeadMin)
 594   {
 595     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
 596     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
 597   }
 598
 599   // start from multiplicity 1
 600   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
 601
 602   for (UInt_t region=0; region<fkRegions; region++)
 603   {
 604     Int_t total = 0;
 605     Int_t count = 0;
 606     Int_t vars[4];
 607
 608     for (Int_t i=0; i<4; i++)
 609       vars[i] = binBegin[i];
 610
 611     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
 612     while (1)
 613     {
 614       if (grid->GetElement(vars) == 0)
 615       {
 616         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
 617           grid->GetBinCenter(0, vars[0]),
 618           grid->GetBinCenter(1, vars[1]),
 619           grid->GetBinCenter(2, vars[2]),
 620           grid->GetBinCenter(3, vars[3])
 621         );
 622         count++;
 623       }
 624
 625       vars[3]++;
 626       for (Int_t i=3; i>0; i--)
 627       {
 628         if (vars[i] == binEnd[i]+1)
 629         {
 630           vars[i] = binBegin[i];
 631           vars[i-1]++;
 632         }
 633       }
 634
 635       if (vars[0] == binEnd[0]+1)
 636         break;
 637       total++;
 638     }
 639
 640     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
 641   }
 642 }
 643
 644 //____________________________________________________________________
 645 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)
 646 {
 647   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 648   //
 649   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 650   // Histogram has to be deleted by the caller of the function
 651   //
 652   // twoD: 0: 1D histogram as function of phi
 653   //       1: 2D histogram, deltaphi, deltaeta
 654   //       10: 1D histogram, within |deltaeta| < 0.8
 655   //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
 656   //
 657   // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
 658   //
 659   // normEvents: if non-0 the number of events/trigger particles for the normalization is filled
 660
 661   // unzoom all axes
 662   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 663   ResetBinLimits(fEventHist->GetGrid(step));
 664
 665   SetBinLimits(fTrackHist[region]->GetGrid(step));
 666
 667   // z vtx
 668   if (fZVtxMax > fZVtxMin)
 669   {
 670     Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
 671     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
 672     fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
 673   }
 674
 675   TH1D* tracks = 0;
 676
 677   if (ptLeadMin < 0)
 678   {
 679     tracks = fTrackHist[region]->ShowProjection(2, step);
 680     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 681     if (fCombineMinMax && region == kMin)
 682     {
 683       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 684       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 685
 686       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 687       tracks->Add(tracks2);
 688
 689       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 690     }
 691
 692     // normalize to get a density (deta dphi)
 693     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 694     Float_t phiRegion = TMath::TwoPi() / 3;
 695     if (!fCombineMinMax && region == kMin)
 696       phiRegion /= 2;
 697     Float_t normalization = phiRegion;
 698     if (etaNorm)
 699       normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 700     //Printf("Normalization: %f", normalization);
 701     tracks->Scale(1.0 / normalization);
 702
 703     TH1D* events = fEventHist->ShowProjection(0, step);
 704     tracks->Divide(events);
 705   }
 706   else
 707   {
 708     if (multBinEnd >= multBinBegin)
 709     {
 710       Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 711       fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 712       fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 713     }
 714
 715     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 716     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 717
 718     Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 719
 720     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
 721
 722     if (twoD == 0)
 723       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
 724     else
 725       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
 726
 727     Printf("Calculated histogram --> %f tracks", tracks->Integral());
 728     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 729
 730     if (twoD == 10 || twoD == 11)
 731     {
 732       Float_t etaLimit = 1.0;
 733       if (twoD == 10)
 734       {
 735         tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
 736
 737         // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
 738         tracks->Scale(1. / 0.75);
 739       }
 740       else if (twoD == 11)
 741       {
 742         TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
 743         TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
 744
 745         tracksTmp1->Add(tracksTmp2);
 746
 747         delete tracks;
 748         tracks = tracksTmp1;
 749         delete tracksTmp2;
 750
 751         tracks->Scale(1. / 0.25);
 752       }
 753     }
 754
 755     // normalize to get a density (deta dphi)
 756     // TODO normalization may be off for 2d histogram
 757     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
 758
 759     if (etaNorm)
 760     {
 761       TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 762       if (strcmp(axis->GetTitle(), "#eta") == 0)
 763       {
 764         Printf("Normalizing using eta-axis range");
 765         normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
 766       }
 767       else
 768       {
 769         Printf("Normalizing assuming accepted range of |eta| < 0.8");
 770         normalization *= 0.8 * 2;
 771       }
 772     }
 773
 774     //Printf("Normalization: %f", normalization);
 775     tracks->Scale(1.0 / normalization);
 776
 777     // 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!
 778     TH1D* events = fEventHist->ShowProjection(0, step);
 779     Int_t nEvents = (Int_t) events->Integral(firstBin, lastBin);
 780     Printf("Calculated histogram --> %d events", nEvents);
 781     if (normEvents)
 782       *normEvents = nEvents;
 783
 784     if (nEvents > 0)
 785       tracks->Scale(1.0 / nEvents);
 786
 787     delete events;
 788   }
 789
 790   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 791   ResetBinLimits(fEventHist->GetGrid(step));
 792
 793   return tracks;
 794 }
 795
 796 //____________________________________________________________________
 797 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)
 798 {
 799   // Calculates a 3d histogram with deltaphi, deltaeta, zvtx on track level and a 1d histogram on event level (as fct of zvtx)
 800   // Histogram has to be deleted by the caller of the function
 801
 802   if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
 803     AliFatal("Histogram without vertex axis provided");
 804
 805   // unzoom all axes
 806   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 807   ResetBinLimits(fEventHist->GetGrid(step));
 808
 809   SetBinLimits(fTrackHist[region]->GetGrid(step));
 810
 811   if (multBinEnd >= multBinBegin)
 812   {
 813     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 814     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 815     fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 816   }
 817
 818   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 819   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 820   Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 821   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
 822   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
 823
 824   *trackHist = (TH3*) fTrackHist[region]->GetGrid(step)->Project(4, 0, 5);
 825   *eventHist = (TH1*) fEventHist->GetGrid(step)->Project(2);
 826
 827   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 828   ResetBinLimits(fEventHist->GetGrid(step));
 829 }
 830
 831 //____________________________________________________________________
 832 void AliUEHist::GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist)
 833 {
 834   // Calculates a 4d histogram with deltaphi, deltaeta, zvtx, multiplicity on track level and
 835   // a 2d histogram on event level (as fct of zvtx, multiplicity)
 836   // Histograms has to be deleted by the caller of the function
 837
 838   if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
 839     AliFatal("Histogram without vertex axis provided");
 840
 841   THnBase* sparse = fTrackHist[region]->GetGrid(step)->GetGrid();
 842   if (fGetMultCacheOn)
 843   {
 844     if (!fGetMultCache)
 845     {
 846       fGetMultCache = ChangeToThn(sparse);
 847       // should work but causes SEGV in ProjectionND below
 848 //       delete sparse;
 849     }
 850     sparse = fGetMultCache;
 851   }
 852
 853   // unzoom all axes
 854   ResetBinLimits(sparse);
 855   ResetBinLimits(fEventHist->GetGrid(step));
 856
 857   SetBinLimits(sparse);
 858
 859   Int_t firstBin = sparse->GetAxis(2)->FindBin(ptLeadMin);
 860   Int_t lastBin = sparse->GetAxis(2)->FindBin(ptLeadMax);
 861   Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 862   sparse->GetAxis(2)->SetRange(firstBin, lastBin);
 863   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
 864
 865   Int_t dimensions[] = { 4, 0, 5, 3 };
 866   THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E");
 867   *eventHist = (TH2*) fEventHist->GetGrid(step)->Project(2, 1);
 868
 869   ResetBinLimits(sparse);
 870   ResetBinLimits(fEventHist->GetGrid(step));
 871
 872   // convert to THn
 873   *trackHist = ChangeToThn(tmpTrackHist);
 874 }
 875
 876 //____________________________________________________________________
 877 TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd)
 878 {
 879   // Calls GetUEHist(...) for *each* vertex bin and multiplicity bin and performs a sum of ratios:
 880   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
 881   // where N is the total number of events/trigger particles and the subscript is the vertex/multiplicity bin
 882   // where mixed is normalized such that the information about the number of pairs in same is kept
 883   //
 884   // returns a 2D histogram: deltaphi, deltaeta
 885   //
 886   // Parameters:
 887   //   mixed: AliUEHist containing mixed event corresponding to this object
 888   //   <other parameters> : check documentation of AliUEHist::GetUEHist
 889
 890   // do not add this hists to the directory
 891   Bool_t oldStatus = TH1::AddDirectoryStatus();
 892   TH1::AddDirectory(kFALSE);
 893
 894   TH2* totalTracks = 0;
 895
 896   THnBase* trackSameAll = 0;
 897   THnBase* trackMixedAll = 0;
 898   TH2* eventSameAll = 0;
 899   TH2* eventMixedAll = 0;
 900
 901   Int_t totalEvents = 0;
 902   Int_t nCorrelationFunctions = 0;
 903
 904   GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackSameAll, &eventSameAll);
 905   mixed->GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackMixedAll, &eventMixedAll);
 906
 907 //   TH1* normParameters = new TH1F("normParameters", "", 100, 0, 2);
 908
 909 //   trackSameAll->Dump();
 910
 911   TAxis* multAxis = trackSameAll->GetAxis(3);
 912
 913   if (multBinEnd < multBinBegin)
 914   {
 915     multBinBegin = 1;
 916     multBinEnd = multAxis->GetNbins();
 917   }
 918
 919   for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++)
 920   {
 921     trackSameAll->GetAxis(3)->SetRange(multBin, multBin);
 922     trackMixedAll->GetAxis(3)->SetRange(multBin, multBin);
 923
 924     // get mixed normalization correction factor: is independent of vertex bin if scaled with number of triggers
 925     trackMixedAll->GetAxis(2)->SetRange(0, -1);
 926     TH2* tracksMixed = trackMixedAll->Projection(1, 0, "E");
 927     Float_t binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1);
 928
 929     // get mixed event normalization by assuming full acceptance at deta at 0 (integrate over dphi), excluding (0, 0)
 930     Double_t mixedNormError = 0;
 931     Double_t mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetYaxis()->FindBin(-0.01)-1, tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
 932     Double_t mixedNormError2 = 0;
 933     Double_t mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetYaxis()->FindBin(0.01)+1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2);
 934
 935     if (mixedNormError == 0 || mixedNormError2 == 0)
 936     {
 937       Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
 938       continue;
 939     }
 940
 941     Int_t nBinsMixedNorm = (tracksMixed->GetYaxis()->FindBin(-0.01) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
 942     mixedNorm /= nBinsMixedNorm;
 943     mixedNormError /= nBinsMixedNorm;
 944
 945     Int_t nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetYaxis()->FindBin(0.01) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
 946     mixedNorm2 /= nBinsMixedNorm2;
 947     mixedNormError2 /= nBinsMixedNorm2;
 948
 949     mixedNorm = mixedNorm / mixedNormError / mixedNormError + mixedNorm2 / mixedNormError2 / mixedNormError2;
 950     mixedNormError = TMath::Sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2));
 951     mixedNorm *= mixedNormError * mixedNormError;
 952
 953 /*    Double_t mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
 954     Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);*/
 955
 956     delete tracksMixed;
 957
 958     Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin);
 959 //     Printf("%f +- %f | %f | %f", mixedNorm, mixedNormError, triggers, mixedNorm / triggers);
 960     if (triggers <= 0)
 961     {
 962       Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
 963       continue;
 964     }
 965
 966     mixedNorm /= triggers;
 967     mixedNormError /= triggers;
 968
 969     if (mixedNorm <= 0)
 970     {
 971       Printf("ERROR: Skipping multiplicity %d because mixed event is empty at (0,0)", multBin);
 972       continue;
 973     }
 974
 975     // finite bin correction
 976     if (fTrackEtaCut > 0)
 977     {
 978       Double_t finiteBinCorrection = -1.0 / (2*fTrackEtaCut) * binWidthEta / 2 + 1;
 979       Printf("Finite bin correction: %f", finiteBinCorrection);
 980       mixedNorm *= finiteBinCorrection;
 981       mixedNormError *= finiteBinCorrection;
 982     }
 983     else
 984     {
 985       Printf("ERROR: fTrackEtaCut not set. Finite bin correction cannot be applied. Continuing anyway...");
 986     }
 987
 988 //     Printf("Norm: %f +- %f", mixedNorm, mixedNormError);
 989
 990 //     normParameters->Fill(mixedNorm);
 991
 992     TAxis* vertexAxis = trackSameAll->GetAxis(2);
 993     for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
 994 //     for (Int_t vertexBin = 3; vertexBin <= 8; vertexBin++)
 995     {
 996       trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
 997       trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
 998
 999       TH2* tracksSame = trackSameAll->Projection(1, 0, "E");
1000       tracksMixed = trackMixedAll->Projection(1, 0, "E");
1001
1002       // asssume flat in dphi, gain in statistics
1003       //     TH1* histMixedproj = mixedTwoD->ProjectionY();
1004       //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
1005       //
1006       //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
1007       //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
1008       //        mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1009
1010       //       delete histMixedproj;
1011
1012       Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin);
1013       if (triggers2 <= 0)
1014       {
1015         Printf("ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin);
1016       }
1017       else
1018       {
1019         tracksMixed->Scale(1.0 / triggers2 / mixedNorm);
1020         // tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
1021
1022 //      new TCanvas; tracksSame->DrawClone("SURF1");
1023 //      new TCanvas; tracksMixed->DrawClone("SURF1");
1024
1025         // some code to judge the relative contribution of the different correlation functions to the overall uncertainty
1026         Double_t sums[] = { 0, 0, 0 };
1027         Double_t errors[] = { 0, 0, 0 };
1028
1029         for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1030           for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1031           {
1032             sums[0] += tracksSame->GetBinContent(x, y);
1033             errors[0] += tracksSame->GetBinError(x, y);
1034             sums[1] += tracksMixed->GetBinContent(x, y);
1035             errors[1] += tracksMixed->GetBinError(x, y);
1036           }
1037
1038         tracksSame->Divide(tracksMixed);
1039
1040         for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1041           for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1042           {
1043             sums[2] += tracksSame->GetBinContent(x, y);
1044             errors[2] += tracksSame->GetBinError(x, y);
1045           }
1046
1047         for (Int_t x=0; x<3; x++)
1048           if (sums[x] > 0)
1049             errors[x] /= sums[x];
1050
1051         Printf("The correlation function %d %d has uncertainties %f %f %f (Ratio S/M %f)", multBin, vertexBin, errors[0], errors[1], errors[2], (errors[1] > 0) ? errors[0] / errors[1] : -1);
1052         // code to draw contributions
1053         /*
1054         TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
1055         proj->SetTitle(Form("Bin %d", vertexBin));
1056         proj->SetLineColor(vertexBin);
1057         proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
1058         */
1059
1060         if (!totalTracks)
1061           totalTracks = (TH2*) tracksSame->Clone("totalTracks");
1062         else
1063           totalTracks->Add(tracksSame);
1064
1065         totalEvents += eventSameAll->GetBinContent(vertexBin, multBin);
1066
1067 //      new TCanvas; tracksMixed->DrawCopy("SURF1");
1068       }
1069
1070       delete tracksSame;
1071       delete tracksMixed;
1072
1073       nCorrelationFunctions++;
1074     }
1075   }
1076
1077   if (totalTracks) {
1078     Double_t sums[] = { 0, 0, 0 };
1079     Double_t errors[] = { 0, 0, 0 };
1080
1081     for (Int_t x=1; x<=totalTracks->GetNbinsX(); x++)
1082       for (Int_t y=1; y<=totalTracks->GetNbinsY(); y++)
1083       {
1084         sums[0] += totalTracks->GetBinContent(x, y);
1085         errors[0] += totalTracks->GetBinError(x, y);
1086       }
1087     if (sums[0] > 0)
1088       errors[0] /= sums[0];
1089
1090     Printf("Dividing %f tracks by %d events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]);
1091     if (totalEvents > 0)
1092       totalTracks->Scale(1.0 / totalEvents);
1093
1094     // normalizate to dphi width
1095     Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1);
1096     totalTracks->Scale(1.0 / normalization);
1097   }
1098
1099   delete trackSameAll;
1100   delete trackMixedAll;
1101   delete eventSameAll;
1102   delete eventMixedAll;
1103
1104 //   new TCanvas; normParameters->Draw();
1105
1106   TH1::AddDirectory(oldStatus);
1107
1108   return totalTracks;
1109 }
1110
1111 /*
1112 TH2* AliUEHist::GetPtDistInPhiRegion(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiBegin, Float_t phiEnd)
1113 {
1114   // Returns pT,assoc distribution in the given pt,trig, multiplicity, phi region
1115   // Does not use sum of ratios for mixed event correction (TODO to be improved)
1116   // returns a 2D histogram: deltaphi, deltaeta
1117   //
1118   // Parameters:
1119   //   mixed: AliUEHist containing mixed event corresponding to this object
1120 */
1121
1122 //____________________________________________________________________
1123 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)
1124 {
1125   // Calls GetUEHist(...) for *each* multiplicity bin and performs a sum of ratios:
1126   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
1127   // where N is the total number of events/trigger particles and the subscript is the multiplicity bin
1128   //
1129   // Can only be used for the 2D histogram at present
1130   //
1131   // Parameters:
1132   //   mixed: AliUEHist containing mixed event corresponding to this object
1133   //   <other parameters> : check documentation of AliUEHist::GetUEHist
1134
1135   TH2* totalTracks = 0;
1136   Int_t totalEvents = 0;
1137
1138   Int_t vertexBin = 1;
1139   TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
1140   if (useVertexBins && !vertexAxis)
1141   {
1142     Printf("Vertex axis requested but not available");
1143     return 0;
1144   }
1145
1146   // vertex bin loop
1147   while (1)
1148   {
1149     if (useVertexBins)
1150     {
1151       SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1152       mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1153       vertexBin++;
1154     }
1155
1156     // multiplicity loop
1157     Int_t multIter = multBinBegin;
1158     while (1)
1159     {
1160       Int_t multBinBeginLocal = multBinBegin;
1161       Int_t multBinEndLocal = multBinEnd;
1162
1163       if (multBinEnd >= multBinBegin)
1164       {
1165         multBinBeginLocal = multIter;
1166         multBinEndLocal = multIter;
1167         multIter++;
1168       }
1169
1170       Int_t nEvents = 0;
1171       TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
1172       // undo normalization
1173       tracks->Scale(nEvents);
1174       totalEvents += nEvents;
1175
1176       TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
1177
1178       // asssume flat in dphi, gain in statistics
1179   //     TH1* histMixedproj = mixedTwoD->ProjectionY();
1180   //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
1181   //
1182   //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
1183   //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
1184   //    mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1185
1186   //       delete histMixedproj;
1187
1188       // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
1189 /*      Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
1190       mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
1191       tracks->Scale(mixedNorm);*/
1192
1193       tracks->Scale(mixedTwoD->Integral() / tracks->Integral());
1194
1195       tracks->Divide(mixedTwoD);
1196
1197       delete mixedTwoD;
1198
1199       if (!totalTracks)
1200         totalTracks = tracks;
1201       else
1202       {
1203         totalTracks->Add(tracks);
1204         delete tracks;
1205       }
1206
1207       if (multIter > multBinEnd)
1208         break;
1209     }
1210
1211     if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
1212       break;
1213   }
1214
1215   if (useVertexBins)
1216     totalEvents = vertexAxis->GetNbins();
1217
1218   Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
1219   if (totalEvents > 0)
1220     totalTracks->Scale(1.0 / totalEvents);
1221
1222   return totalTracks;
1223 }
1224
1225 //____________________________________________________________________
1226 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)
1227 {
1228   // returns a histogram projected to pT,assoc with the provided cuts
1229
1230    // unzoom all axes
1231   ResetBinLimits(fTrackHist[region]->GetGrid(step));
1232   ResetBinLimits(fEventHist->GetGrid(step));
1233
1234   TH1D* tracks = 0;
1235
1236     // 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
1237   // therefore the number density must be calculated here per leading pT bin and then summed
1238
1239   if (multBinEnd > multBinBegin)
1240     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
1241   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
1242   fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
1243
1244   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
1245   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
1246
1247   // get real phi cuts due to binning
1248   Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
1249   Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
1250   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);
1251
1252   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
1253   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
1254
1255   TH1D* events = fEventHist->ShowProjection(0, step);
1256
1257   for (Int_t bin=firstBin; bin<=lastBin; bin++)
1258   {
1259     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
1260
1261     // project to pT,assoc
1262     TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
1263
1264     Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
1265     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
1266
1267     // normalize to get a density (deta dphi)
1268     Float_t normalization = 1;
1269     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
1270     if (strcmp(axis->GetTitle(), "#eta") == 0)
1271     {
1272       Printf("Normalizing using eta-axis range");
1273       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
1274     }
1275     else
1276     {
1277       Printf("Normalizating assuming accepted range of |eta| < 0.8");
1278       normalization *= 0.8 * 2;
1279     }
1280
1281     // dphi
1282     if (!skipPhiNormalization)
1283       normalization *= phiMaxReal - phiMinReal;
1284
1285     //Printf("Normalization: %f", normalization);
1286     tracksTmp->Scale(1.0 / normalization);
1287
1288     // and now dpT (bins have different width)
1289     for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
1290     {
1291       tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1292       tracksTmp->SetBinError  (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1293     }
1294
1295     Int_t nEvents = (Int_t) events->GetBinContent(bin);
1296     if (nEvents > 0)
1297       tracksTmp->Scale(1.0 / nEvents);
1298     Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
1299
1300     if (!tracks)
1301       tracks = tracksTmp;
1302     else
1303     {
1304       tracks->Add(tracksTmp);
1305       delete tracksTmp;
1306     }
1307   }
1308
1309   delete events;
1310
1311   ResetBinLimits(fTrackHist[region]->GetGrid(step));
1312   ResetBinLimits(fEventHist->GetGrid(step));
1313
1314   return tracks;
1315 }
1316
1317 void AliUEHist::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
1318 {
1319   // multiplies <grid> with <histogram> and put the result in <target>
1320   // <grid> has usually more dimensions than histogram. The axis which are used to choose the value
1321   // from <histogram> are given in <var1> and <var2>
1322   //
1323   // if <histogram> is 0, just copies content from step1 to step2
1324
1325   // clear target histogram
1326   if (grid != target)
1327     target->Reset();
1328
1329   if (histogram != 0)
1330   {
1331     if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
1332       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
1333
1334     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
1335       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
1336   }
1337
1338   if (grid->GetNdimensions() > 6)
1339     AliFatal("Too many dimensions in THnSparse");
1340
1341   Int_t bins[6];
1342
1343   // optimized implementation
1344   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1345   {
1346     Double_t value = grid->GetBinContent(binIdx, bins);
1347     Double_t error = grid->GetBinError(binIdx);
1348
1349     if (histogram != 0)
1350     {
1351       if (var2 < 0)
1352       {
1353         value *= histogram->GetBinContent(bins[var1]);
1354         error *= histogram->GetBinContent(bins[var1]);
1355       }
1356       else
1357       {
1358         value *= histogram->GetBinContent(bins[var1], bins[var2]);
1359         error *= histogram->GetBinContent(bins[var1], bins[var2]);
1360       }
1361     }
1362
1363     target->SetBinContent(bins, value);
1364     target->SetBinError(bins, error);
1365   }
1366 }
1367
1368 //____________________________________________________________________
1369 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
1370 {
1371   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
1372   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
1373   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
1374   //
1375   // if trackCorrection is 0, just copies content from step1 to step2
1376
1377   for (UInt_t region=0; region<fkRegions; region++)
1378     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
1379 }
1380
1381 //____________________________________________________________________
1382 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
1383 {
1384   //
1385   // see documentation of CorrectTracks above
1386   //
1387
1388   if (!fTrackHist[region])
1389     return;
1390
1391   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
1392   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
1393
1394   Float_t entriesBefore = grid->GetEntries();
1395
1396   MultiplyHistograms(grid, target, trackCorrection, var1, var2);
1397
1398   Printf("AliUEHist::CorrectTracks: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0);
1399 }
1400
1401 //____________________________________________________________________
1402 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
1403 {
1404   // corrects from step1 to step2 by multiplying the events with eventCorrection
1405   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
1406   //
1407   // if eventCorrection is 0, just copies content from step1 to step2
1408
1409   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
1410   AliCFGridSparse* target = fEventHist->GetGrid(step2);
1411
1412   Float_t entriesBefore = grid->GetEntries();
1413
1414   MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);
1415
1416   Printf("AliUEHist::CorrectEvents: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0);
1417 }
1418
1419 //____________________________________________________________________
1420 void AliUEHist::Correct(AliUEHist* corrections)
1421 {
1422   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
1423   //
1424   // in this object the data is expected in the step kCFStepReconstructed
1425
1426   if (fHistogramType.Length() == 0)
1427   {
1428     Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
1429     if (fTrackHist[kToward]->GetNVar() < 5)
1430     {
1431       if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#varphid#eta") == 0)
1432         fHistogramType = "NumberDensitypT";
1433       else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#eta") == 0)
1434         fHistogramType = "SumpT";
1435     }
1436     else if (fTrackHist[kToward]->GetNVar() == 5)
1437     {
1438       if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
1439         fHistogramType = "NumberDensityPhi";
1440       else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
1441         fHistogramType = "NumberDensityPhiCentrality";
1442     }
1443
1444     if (fHistogramType.Length() == 0)
1445       AliFatal("Cannot figure out histogram type. Try setting it manually...");
1446   }
1447
1448   Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
1449
1450   if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
1451   {
1452     // ---- track level
1453
1454     // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
1455     // extracted as function of leading pT
1456     Bool_t biasFromMC = kFALSE;
1457     const char* projAxis = "z";
1458     Int_t secondBin = -1;
1459
1460     if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
1461     {
1462       projAxis = "yz";
1463       secondBin = 4;
1464     }
1465
1466     for (UInt_t region = 0; region < fkRegions; region++)
1467     {
1468       if (!fTrackHist[region])
1469         continue;
1470
1471       TH1* leadingBiasTracks =  0;
1472       if (biasFromMC)
1473       {
1474         leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
1475         Printf("WARNING: Using MC bias correction");
1476       }
1477       else
1478         leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
1479
1480       CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
1481       if (region == kMin && fCombineMinMax)
1482       {
1483         CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
1484         delete leadingBiasTracks;
1485         break;
1486       }
1487       delete leadingBiasTracks;
1488     }
1489
1490     TH1* leadingBiasEvents = 0;
1491     if (biasFromMC)
1492       leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
1493     else
1494       leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
1495
1496     //new TCanvas; leadingBiasEvents->DrawCopy();
1497     CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
1498     delete leadingBiasEvents;
1499
1500     // --- contamination correction ---
1501
1502     TH2D* contamination = corrections->GetTrackingContamination();
1503     if (corrections->fContaminationEnhancement)
1504     {
1505       Printf("Applying contamination enhancement");
1506
1507       for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
1508         for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
1509           contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
1510     }
1511     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
1512     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1513     delete contamination;
1514
1515     // --- efficiency correction ---
1516     // correct single-particle efficiency for associated particles
1517     // in addition correct for efficiency on trigger particles (tracks AND events)
1518
1519     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
1520     // use kCFStepVertex as a temporary step (will be overwritten later anyway)
1521     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1522     delete efficiencyCorrection;
1523
1524     // correct pT,T
1525     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
1526     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
1527     CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
1528     delete efficiencyCorrection;
1529
1530     // copy
1531     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
1532     CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
1533
1534     // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
1535     // practically independent of low pT cut
1536     TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
1537
1538     // convert stage from true multiplicity to observed multiplicity by simple conversion factor
1539     TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
1540     vertexCorrectionObs->Reset();
1541
1542     TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
1543     // some defaults
1544     func->SetParameters(0.1, 1, -0.7);
1545     vertexCorrection->Fit(func, "0I", "", 0, 3);
1546     for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
1547     {
1548       Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
1549       if (xPos < 3)
1550         vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
1551       else
1552         vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
1553     }
1554
1555     #if 0
1556     new TCanvas;
1557     vertexCorrection->DrawCopy();
1558     vertexCorrectionObs->SetLineColor(2);
1559     vertexCorrectionObs->DrawCopy("same");
1560     func->SetRange(0, 4);
1561     func->DrawClone("same");
1562     #endif
1563
1564     CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
1565     CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
1566     delete vertexCorrectionObs;
1567     delete vertexCorrection;
1568     delete func;
1569
1570     // copy
1571     CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
1572     CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
1573   }
1574   else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
1575   {
1576     if (fTrackHist[0]->GetNVar() <= 5)
1577     {
1578       // do corrections copying between steps
1579 //       CFStep step = kCFStepReconstructed;
1580       CFStep step = kCFStepBiasStudy;
1581
1582       // copy
1583       CorrectTracks(step, kCFStepTracked, 0, -1);
1584       CorrectEvents(step, kCFStepTracked, 0, -1);
1585
1586       // Dont use eta in the following, because it is a Delta-eta axis
1587
1588       // contamination correction
1589       // correct single-particle contamination for associated particles
1590
1591       TH1* contamination = corrections->GetTrackingContamination(1);
1592
1593       if (0)
1594       {
1595         Printf("Applying contamination enhancement");
1596
1597         for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1598         {
1599           printf("%f", contamination->GetBinContent(bin));
1600           if (contamination->GetBinContent(bin) > 0)
1601             contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1602           printf(" --> %f\n", contamination->GetBinContent(bin));
1603         }
1604       }
1605
1606       CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1607       delete contamination;
1608
1609       // correct for additional contamination due to trigger particle around phi ~ 0
1610       TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1611       if (0)
1612       {
1613         Printf("Applying contamination enhancement");
1614
1615         for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1616           for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1617           {
1618             printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1619             if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1620               correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1621             printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1622           }
1623       }
1624
1625 //       new TCanvas; correlatedContamination->DrawCopy("COLZ");
1626       CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1627 //       Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1628
1629       delete correlatedContamination;
1630
1631       // TODO correct for contamination of trigger particles (for tracks AND events)
1632       CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1633
1634       // --- efficiency correction ---
1635       // correct single-particle efficiency for associated particles
1636       // in addition correct for efficiency on trigger particles (tracks AND events)
1637
1638       // in bins of pT and centrality
1639       TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1640 //       new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1641       // use kCFStepAnaTopology as a temporary step
1642       CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1643       delete efficiencyCorrection;
1644
1645       // correct pT,T in bins of pT and centrality
1646       efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1647       CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1648       CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1649       delete efficiencyCorrection;
1650
1651       // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1652       // copy
1653       CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1654       CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1655     }
1656     else
1657     {
1658       // with 6 axes there is not enough memory, do the corrections in-place
1659       Printf(">>>>>>>> Applying corrections in place to reduce memory consumption");
1660       CFStep step = kCFStepBiasStudy;
1661 //       CFStep step = kCFStepReconstructed;
1662
1663       // Dont use eta in the following, because it is a Delta-eta axis
1664
1665       // --- contamination correction ---
1666       // correct single-particle contamination for associated particles
1667       // correct contamination for trigger particles (tracks AND events)
1668
1669       // in bins of p,T
1670       TH1* contamination = corrections->GetTrackingContamination(1);
1671
1672       if (0)
1673       {
1674         Printf("Applying contamination enhancement");
1675
1676         for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1677         {
1678           printf("%f", contamination->GetBinContent(bin));
1679           if (contamination->GetBinContent(bin) > 0)
1680             contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1681           printf(" --> %f\n", contamination->GetBinContent(bin));
1682         }
1683       }
1684
1685       // correct pT,A in bins of pT
1686       CorrectTracks(step, step, contamination, 1);
1687       delete contamination;
1688
1689       // correct pT,T in bins of pT (for tracks AND events)
1690       contamination = corrections->GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, -1, 2);
1691       new TCanvas; contamination->DrawCopy();
1692       CorrectEvents(step, step, contamination, 0);
1693       CorrectTracks(step, step, contamination, 2);
1694       delete contamination;
1695
1696       // correct for additional contamination due to trigger particle around phi ~ 0
1697       if (0)
1698       {
1699         TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1700         if (1)
1701         {
1702           Printf("Applying contamination enhancement");
1703
1704           for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1705             for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1706             {
1707               printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1708               if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1709                 correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1710               printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1711             }
1712         }
1713
1714         // new TCanvas; correlatedContamination->DrawCopy("COLZ");
1715         CorrectCorrelatedContamination(step, 0, correlatedContamination);
1716
1717         delete correlatedContamination;
1718       }
1719       else
1720         Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1721
1722       // --- tracking efficiency correction ---
1723       // correct single-particle efficiency for associated particles
1724       // correct for efficiency on trigger particles (tracks AND events)
1725
1726       // in bins of pT and centrality
1727       TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1728       // new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1729
1730       // correct pT,A in bins of pT and centrality
1731       CorrectTracks(step, step, efficiencyCorrection, 1, 3);
1732       delete efficiencyCorrection;
1733
1734       // correct pT,T in bins of pT and centrality
1735       efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1736       CorrectEvents(step, step, efficiencyCorrection, 0, 1);
1737       CorrectTracks(step, step, efficiencyCorrection, 2, 3);
1738
1739       delete efficiencyCorrection;
1740     }
1741   }
1742   else
1743     AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1744 }
1745
1746 //____________________________________________________________________
1747 THnBase* AliUEHist::GetTrackEfficiencyND(CFStep step1, CFStep step2)
1748 {
1749   // creates a track-level efficiency by dividing step2 by step1
1750   // in all dimensions but the particle species one
1751
1752   AliCFContainer* sourceContainer = fTrackHistEfficiency;
1753   // step offset because we start with kCFStepAnaTopology
1754   step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1755   step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1756
1757   ResetBinLimits(sourceContainer->GetGrid(step1));
1758   ResetBinLimits(sourceContainer->GetGrid(step2));
1759
1760   if (fEtaMax > fEtaMin)
1761   {
1762     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1763     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1764     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1765   }
1766
1767   Int_t dimensions[] = { 0, 1, 3, 4 };
1768   THnBase* generated = sourceContainer->GetGrid(step1)->GetGrid()->ProjectionND(4, dimensions);
1769   THnBase* measured = sourceContainer->GetGrid(step2)->GetGrid()->ProjectionND(4, dimensions);
1770
1771 //   Printf("%d %d %f %f", step1, step2, generated->GetEntries(), measured->GetEntries());
1772
1773   ResetBinLimits(sourceContainer->GetGrid(step1));
1774   ResetBinLimits(sourceContainer->GetGrid(step2));
1775
1776   measured->Divide(measured, generated, 1, 1, "B");
1777
1778   delete generated;
1779
1780   return measured;
1781 }
1782
1783 //____________________________________________________________________
1784 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1785 {
1786   // creates a track-level efficiency by dividing step2 by step1
1787   // projected to axis1 and axis2 (optional if >= 0)
1788   //
1789   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1790
1791   // integrate over regions
1792   // cache it for efficiency (usually more than one efficiency is requested)
1793
1794   AliCFContainer* sourceContainer = 0;
1795
1796   if (source == 0)
1797   {
1798     if (!fCache)
1799     {
1800       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1801       for (UInt_t i = 1; i < fkRegions; i++)
1802         if (fTrackHist[i])
1803           fCache->Add(fTrackHist[i]);
1804     }
1805     sourceContainer = fCache;
1806   }
1807   else if (source == 1 || source == 2)
1808   {
1809     sourceContainer = fTrackHistEfficiency;
1810     // step offset because we start with kCFStepAnaTopology
1811     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1812     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1813   }
1814   else
1815     return 0;
1816
1817   // reset all limits and set the right ones except those in axis1, axis2 and axis3
1818   ResetBinLimits(sourceContainer->GetGrid(step1));
1819   ResetBinLimits(sourceContainer->GetGrid(step2));
1820   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1821   {
1822     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1823     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1824     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1825   }
1826   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1827   {
1828     Printf("Restricted pt-range to %f %f", fPtMin, fPtMax);
1829     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1830     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1831   }
1832   if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1833   {
1834     Printf("Restricted centrality range to %f %f", fCentralityMin, fCentralityMax);
1835     sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1836     sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1837   }
1838   if (fZVtxMax > fZVtxMin && axis1 != 4 && axis2 != 4 && axis3 != 4)
1839   {
1840     Printf("Restricted z-vtx range to %f %f", fZVtxMin, fZVtxMax);
1841     sourceContainer->GetGrid(step1)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1842     sourceContainer->GetGrid(step2)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1843   }
1844
1845   TH1* measured = 0;
1846   TH1* generated = 0;
1847
1848   if (axis3 >= 0)
1849   {
1850     generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1851     measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1852   }
1853   else if (axis2 >= 0)
1854   {
1855     generated = sourceContainer->Project(step1, axis1, axis2);
1856     measured = sourceContainer->Project(step2, axis1, axis2);
1857   }
1858   else
1859   {
1860     generated = sourceContainer->Project(step1, axis1);
1861     measured = sourceContainer->Project(step2, axis1);
1862   }
1863
1864   // check for bins with less than 50 entries, print warning
1865   Int_t binBegin[3];
1866   Int_t binEnd[3];
1867
1868   binBegin[0] = 1;
1869   binBegin[1] = 1;
1870   binBegin[2] = 1;
1871
1872   binEnd[0] = generated->GetNbinsX();
1873   binEnd[1] = generated->GetNbinsY();
1874   binEnd[2] = generated->GetNbinsZ();
1875
1876   if (fEtaMax > fEtaMin)
1877   {
1878     if (axis1 == 0)
1879     {
1880       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1881       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
1882     }
1883     if (axis2 == 0)
1884     {
1885       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1886       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
1887     }
1888     if (axis3 == 0)
1889     {
1890       binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1891       binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
1892     }
1893   }
1894
1895   if (fPtMax > fPtMin)
1896   {
1897     // TODO this is just checking up to 15 for now
1898     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1899     if (axis1 == 1)
1900     {
1901       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1902       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
1903     }
1904     if (axis2 == 1)
1905     {
1906       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1907       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
1908     }
1909     if (axis3 == 1)
1910     {
1911       binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1912       binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
1913     }
1914   }
1915
1916   Int_t total = 0;
1917   Int_t count = 0;
1918   Int_t vars[3];
1919
1920   for (Int_t i=0; i<3; i++)
1921     vars[i] = binBegin[i];
1922
1923   const Int_t limit = 50;
1924   while (1)
1925   {
1926     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1927     {
1928       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1929       count++;
1930     }
1931     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1932     {
1933       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
1934         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1935         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1936         generated->GetBinContent(vars[0], vars[1]));
1937       count++;
1938     }
1939     else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1940     {
1941       Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
1942         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1943         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1944         generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1945         generated->GetBinContent(vars[0], vars[1], vars[2]));
1946       count++;
1947     }
1948
1949     vars[2]++;
1950     if (vars[2] == binEnd[2]+1)
1951     {
1952       vars[2] = binBegin[2];
1953       vars[1]++;
1954     }
1955
1956     if (vars[1] == binEnd[1]+1)
1957     {
1958       vars[1] = binBegin[1];
1959       vars[0]++;
1960     }
1961
1962     if (vars[0] == binEnd[0]+1)
1963       break;
1964     total++;
1965   }
1966
1967   Printf("Correction has %d empty bins (out of %d bins)", count, total);
1968
1969   // rebin if required
1970   if (source == 2)
1971   {
1972     TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1973
1974     if (axis->GetNbins() < measured->GetNbinsX())
1975     {
1976       if (axis2 != -1)
1977       {
1978         // 2d rebin with variable axis does not exist in root
1979
1980         TH1* tmp = measured;
1981         measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1982         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1983           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1984           {
1985             ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1986             measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1987           }
1988         delete tmp;
1989
1990         tmp = generated;
1991         generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1992         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1993           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1994           {
1995             ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1996             generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1997           }
1998         delete tmp;
1999       }
2000       else
2001       {
2002         TH1* tmp = measured;
2003         measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2004         delete tmp;
2005
2006         tmp = generated;
2007         generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2008         delete tmp;
2009       }
2010     }
2011     else if (axis->GetNbins() > measured->GetNbinsX())
2012     {
2013       if (axis2 != -1)
2014         AliFatal("Rebinning only works for 1d at present");
2015
2016       // this is an unfortunate case where the number of bins has to be increased in principle
2017       // there is a region where the binning is finner in one histogram and a region where it is the other way round
2018       // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
2019       // only a certain binning is supported here
2020       if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
2021         AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
2022
2023       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};
2024
2025       // reduce binning below 5 GeV/c
2026       TH1* tmp = measured;
2027       measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2028       delete tmp;
2029
2030       // increase binning above 5 GeV/c
2031       tmp = measured;
2032       measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2033       for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
2034       {
2035         measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
2036         measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
2037       }
2038       delete tmp;
2039
2040       // reduce binning below 5 GeV/c
2041       tmp = generated;
2042       generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2043       delete tmp;
2044
2045       // increase binning above 5 GeV/c
2046       tmp = generated;
2047       generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2048       for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
2049       {
2050         generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
2051         generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
2052       }
2053       delete tmp;
2054     }
2055   }
2056
2057   measured->Divide(measured, generated, 1, 1, "B");
2058
2059   delete generated;
2060
2061   ResetBinLimits(sourceContainer->GetGrid(step1));
2062   ResetBinLimits(sourceContainer->GetGrid(step2));
2063
2064   return measured;
2065 }
2066
2067 //____________________________________________________________________
2068 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
2069 {
2070   // creates a event-level efficiency by dividing step2 by step1
2071   // projected to axis1 and axis2 (optional if >= 0)
2072
2073   if (ptLeadMax > ptLeadMin)
2074   {
2075     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2076     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2077   }
2078
2079   TH1* measured = 0;
2080   TH1* generated = 0;
2081
2082   if (axis2 >= 0)
2083   {
2084     generated = fEventHist->Project(step1, axis1, axis2);
2085     measured = fEventHist->Project(step2, axis1, axis2);
2086   }
2087   else
2088   {
2089     generated = fEventHist->Project(step1, axis1);
2090     measured = fEventHist->Project(step2, axis1);
2091   }
2092
2093   measured->Divide(measured, generated, 1, 1, "B");
2094
2095   delete generated;
2096
2097   if (ptLeadMax > ptLeadMin)
2098   {
2099     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
2100     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
2101   }
2102
2103   return measured;
2104 }
2105
2106 //____________________________________________________________________
2107 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
2108 {
2109   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
2110   // where the matching is done of the z axis of hist1 with the x axis of hist2
2111
2112   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
2113     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
2114
2115   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
2116   {
2117     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
2118     {
2119       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
2120       {
2121         if (hist2->GetBinContent(z) > 0)
2122         {
2123           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
2124           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
2125         }
2126         else
2127         {
2128           hist1->SetBinContent(x, y, z, 0);
2129           hist1->SetBinError(x, y, z, 0);
2130         }
2131       }
2132     }
2133   }
2134 }
2135
2136 //____________________________________________________________________
2137 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
2138 {
2139   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
2140   // in the given region (sum over all regions is calculated if region == -1)
2141   // done by weighting the track-level distribution with the number of events as function of leading pT
2142   // and then calculating the ratio between the distributions
2143   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
2144   //   no projection is done if axis == 0
2145   // weighting: 0 = tracks weighted with events (as discussed above)
2146   //            1 = only track bias is returned
2147   //            2 = only event bias is returned
2148
2149   AliCFContainer* tmp = 0;
2150
2151   if (region == -1)
2152   {
2153     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
2154     for (UInt_t i = 1; i < fkRegions; i++)
2155       if (fTrackHist[i])
2156         tmp->Add(fTrackHist[i]);
2157   }
2158   else if (region == kMin && fCombineMinMax)
2159   {
2160     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
2161     tmp->Add(fTrackHist[kMax]);
2162   }
2163   else
2164     tmp = fTrackHist[region];
2165
2166   ResetBinLimits(tmp->GetGrid(step1));
2167   ResetBinLimits(fEventHist->GetGrid(step1));
2168   SetBinLimits(tmp->GetGrid(step1));
2169
2170   ResetBinLimits(tmp->GetGrid(step2));
2171   ResetBinLimits(fEventHist->GetGrid(step2));
2172   SetBinLimits(tmp->GetGrid(step2));
2173
2174   TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
2175   TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
2176   if (weighting == 0)
2177     WeightHistogram(hist1, events1);
2178
2179   TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
2180   TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
2181   if (weighting == 0)
2182     WeightHistogram(hist2, events2);
2183
2184   TH1* generated = hist1;
2185   TH1* measured = hist2;
2186
2187   if (weighting == 0 || weighting == 1)
2188   {
2189     if (axis)
2190     {
2191       if (leadPtMax > leadPtMin)
2192       {
2193         hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2194         hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2195       }
2196
2197       if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
2198       {
2199         hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2200         hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2201       }
2202
2203       generated = hist1->Project3D(axis);
2204       measured  = hist2->Project3D(axis);
2205
2206       // delete hists here if projection has been done
2207       delete hist1;
2208       delete hist2;
2209     }
2210     delete events1;
2211     delete events2;
2212   }
2213   else if (weighting == 2)
2214   {
2215     delete hist1;
2216     delete hist2;
2217     generated = events1;
2218     measured = events2;
2219   }
2220
2221   measured->Divide(generated);
2222
2223   delete generated;
2224
2225   ResetBinLimits(tmp->GetGrid(step1));
2226   ResetBinLimits(tmp->GetGrid(step2));
2227
2228   if ((region == -1) || (region == kMin && fCombineMinMax))
2229     delete tmp;
2230
2231   return measured;
2232 }
2233
2234 //____________________________________________________________________
2235 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
2236 {
2237   // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
2238
2239   if (!fTrackHist[region])
2240     return;
2241
2242   THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
2243
2244   Int_t var1 = 1;
2245   Int_t var2 = 2;
2246
2247   if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
2248     AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
2249
2250   if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
2251     AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
2252
2253   // optimized implementation
2254   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2255   {
2256     Int_t bins[6];
2257
2258     Double_t value = grid->GetBinContent(binIdx, bins);
2259     Double_t error = grid->GetBinError(binIdx);
2260
2261     // check eta and phi axes
2262     if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
2263       continue;
2264     if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
2265       continue;
2266
2267     value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2268     error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2269
2270     grid->SetBinContent(bins, value);
2271     grid->SetBinError(bins, error);
2272   }
2273
2274   Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
2275 }
2276
2277 //____________________________________________________________________
2278 TH2* AliUEHist::GetCorrelatedContamination()
2279 {
2280   // 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!)
2281
2282   Int_t step1 = kCFStepTrackedOnlyPrim;
2283   Int_t step2 = kCFStepTracked;
2284
2285   fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
2286   fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
2287   TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
2288
2289   fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
2290   fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
2291   TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
2292
2293   tracksStep1->Divide(tracksStep2);
2294
2295   TH1* triggersStep1 = fEventHist->Project(step1, 0);
2296   TH1* triggersStep2 = fEventHist->Project(step2, 0);
2297
2298   TH1* singleParticle = GetTrackingContamination(1);
2299
2300   for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
2301     for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
2302       if (singleParticle->GetBinContent(x) > 0 && triggersStep1->GetBinContent(y) > 0)
2303         tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
2304       else
2305         tracksStep1->SetBinContent(x, y, 0);
2306
2307   delete singleParticle;
2308   delete tracksStep2;
2309   delete triggersStep1;
2310   delete triggersStep2;
2311
2312   return tracksStep1;
2313 }
2314
2315 //____________________________________________________________________
2316 TH2D* AliUEHist::GetTrackingEfficiency()
2317 {
2318   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2319   // integrates over the regions and all other variables than pT and eta to increase the statistics
2320   //
2321   // returned histogram has to be deleted by the user
2322
2323   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
2324 }
2325
2326 //____________________________________________________________________
2327 TH2D* AliUEHist::GetFakeRate()
2328 {
2329   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepReconstructed, 0, 1));
2330 }
2331
2332 //____________________________________________________________________
2333 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
2334 {
2335   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2336   // integrates over the regions and all other variables than pT, centrality to increase the statistics
2337   //
2338   // returned histogram has to be deleted by the user
2339
2340   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
2341 }
2342
2343 //____________________________________________________________________
2344 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
2345 {
2346   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2347   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2348
2349   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
2350 }
2351
2352 //____________________________________________________________________
2353 TH1D* AliUEHist::GetFakeRate(Int_t axis)
2354 {
2355   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepReconstructed, axis));
2356 }
2357 //____________________________________________________________________
2358 TH2D* AliUEHist::GetTrackingCorrection()
2359 {
2360   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2361   // integrates over the regions and all other variables than pT and eta to increase the statistics
2362   //
2363   // returned histogram has to be deleted by the user
2364
2365   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
2366 }
2367
2368 //____________________________________________________________________
2369 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
2370 {
2371   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2372   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2373
2374   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
2375 }
2376
2377 //____________________________________________________________________
2378 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
2379 {
2380   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2381   // integrates over the regions and all other variables than pT and eta to increase the statistics
2382   //
2383   // returned histogram has to be deleted by the user
2384
2385   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
2386 }
2387
2388 //____________________________________________________________________
2389 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
2390 {
2391   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2392   // integrates over the regions and all other variables than pT and centrality to increase the statistics
2393   //
2394   // returned histogram has to be deleted by the user
2395
2396   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
2397 }
2398
2399 //____________________________________________________________________
2400 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
2401 {
2402   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2403   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2404
2405   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
2406 }
2407
2408 //____________________________________________________________________
2409 TH2D* AliUEHist::GetTrackingContamination()
2410 {
2411   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2412   // integrates over the regions and all other variables than pT and eta to increase the statistics
2413   //
2414   // returned histogram has to be deleted by the user
2415
2416   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
2417 }
2418
2419 //____________________________________________________________________
2420 TH2D* AliUEHist::GetTrackingContaminationCentrality()
2421 {
2422   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2423   // integrates over the regions and all other variables than pT and centrality to increase the statistics
2424   //
2425   // returned histogram has to be deleted by the user
2426
2427   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
2428 }
2429
2430 //____________________________________________________________________
2431 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
2432 {
2433   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2434   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2435
2436   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
2437 }
2438
2439 //____________________________________________________________________
2440 const char* AliUEHist::GetRegionTitle(Region region)
2441 {
2442   // returns the name of the given region
2443
2444   switch (region)
2445   {
2446     case kToward:
2447       return "Towards";
2448     case kAway:
2449       return "Away";
2450     case kMin:
2451       return (fCombineMinMax) ? "Transverse" : "Min";
2452     case kMax:
2453       return "Max";
2454   }
2455
2456   return 0;
2457 }
2458
2459 //____________________________________________________________________
2460 const char* AliUEHist::GetStepTitle(CFStep step)
2461 {
2462   // returns the name of the given step
2463
2464   switch (step)
2465   {
2466     case kCFStepAll:
2467       return "All events";
2468     case kCFStepTriggered:
2469       return "Triggered";
2470     case kCFStepVertex:
2471       return "Primary Vertex";
2472     case kCFStepAnaTopology:
2473       return "Required analysis topology";
2474     case kCFStepTrackedOnlyPrim:
2475       return "Tracked (matched MC, only primaries)";
2476     case kCFStepTracked:
2477       return "Tracked (matched MC, all)";
2478     case kCFStepReconstructed:
2479       return "Reconstructed";
2480     case kCFStepRealLeading:
2481       return "Correct leading particle identified";
2482     case kCFStepBiasStudy:
2483       return "Bias study applying tracking efficiency";
2484     case kCFStepBiasStudy2:
2485       return "Bias study applying tracking efficiency in two steps";
2486     case kCFStepCorrected:
2487       return "Corrected for efficiency on-the-fly";
2488   }
2489
2490   return 0;
2491 }
2492
2493 //____________________________________________________________________
2494 void AliUEHist::CopyReconstructedData(AliUEHist* from)
2495 {
2496   // copies those histograms extracted from ESD to this object
2497
2498   // TODO at present only the pointers are copied
2499
2500   for (Int_t region=0; region<4; region++)
2501   {
2502     if (!fTrackHist[region])
2503       continue;
2504
2505     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
2506     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2507     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
2508   }
2509
2510   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
2511   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2512   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
2513 }
2514
2515 void AliUEHist::DeepCopy(AliUEHist* from)
2516 {
2517   // copies the entries of this object's members from the object <from> to this object
2518   // fills using the fill function and thus allows that the objects have different binning
2519
2520   for (Int_t region=0; region<4; region++)
2521   {
2522     if (!fTrackHist[region] || !from->fTrackHist[region])
2523       continue;
2524
2525     for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
2526     {
2527       Printf("Copying region %d step %d", region, step);
2528       THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
2529       THnSparse* source = from->fTrackHist[region]->GetGrid(step)->GetGrid();
2530
2531       target->Reset();
2532       target->RebinnedAdd(source);
2533     }
2534   }
2535
2536   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2537   {
2538     Printf("Ev: Copying step %d", step);
2539     THnSparse* target = fEventHist->GetGrid(step)->GetGrid();
2540     THnSparse* source = from->fEventHist->GetGrid(step)->GetGrid();
2541
2542     target->Reset();
2543     target->RebinnedAdd(source);
2544   }
2545
2546   for (Int_t step=0; step<TMath::Min(fTrackHistEfficiency->GetNStep(), from->fTrackHistEfficiency->GetNStep()); step++)
2547   {
2548     if (!from->fTrackHistEfficiency->GetGrid(step))
2549       continue;
2550
2551     Printf("Eff: Copying step %d", step);
2552     THnSparse* target = fTrackHistEfficiency->GetGrid(step)->GetGrid();
2553     THnSparse* source = from->fTrackHistEfficiency->GetGrid(step)->GetGrid();
2554
2555     target->Reset();
2556     target->RebinnedAdd(source);
2557   }
2558 }
2559
2560 //____________________________________________________________________
2561 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
2562 {
2563   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
2564
2565   Float_t fitRangeBegin = 5.01;
2566   Float_t fitRangeEnd = 14.99;
2567   Float_t extendRangeBegin = 10.01;
2568
2569   if (fTrackHistEfficiency->GetNVar() == 3)
2570   {
2571     TH1* obj = GetTrackingEfficiency(1);
2572
2573     if (verbose)
2574     {
2575       new TCanvas;
2576       obj->Draw();
2577     }
2578
2579     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2580
2581     Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
2582
2583     obj = GetTrackingContamination(1);
2584
2585     if (verbose)
2586     {
2587       new TCanvas;
2588       obj->Draw();
2589     }
2590
2591     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2592
2593     Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
2594
2595     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);
2596
2597     // extend for full pT range
2598     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2599       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2600         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2601         {
2602
2603           Int_t bins[3];
2604           bins[0] = x;
2605           bins[1] = y;
2606           bins[2] = z;
2607
2608           fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2609           fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
2610           fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
2611         }
2612   }
2613   else if (fTrackHistEfficiency->GetNVar() == 4)
2614   {
2615     // fit in centrality intervals of 20% for efficiency, one bin for contamination
2616     Float_t* trackingEff = 0;
2617     Float_t* trackingCont = 0;
2618     Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
2619     Int_t nCentralityBins = 5;
2620
2621     Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
2622
2623     // 0 = eff; 1 = cont
2624     for (Int_t caseNo = 0; caseNo < 2; caseNo++)
2625     {
2626       Float_t* target = 0;
2627       Int_t centralityBinsLocal = nCentralityBins;
2628
2629       if (caseNo == 0)
2630       {
2631         trackingEff = new Float_t[centralityBinsLocal];
2632         target = trackingEff;
2633       }
2634       else
2635       {
2636         centralityBinsLocal = 1;
2637         trackingCont = new Float_t[centralityBinsLocal];
2638         target = trackingCont;
2639       }
2640
2641       for (Int_t i=0; i<centralityBinsLocal; i++)
2642       {
2643         if (centralityBinsLocal == 1)
2644           SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
2645         else
2646           SetCentralityRange(centralityBins[i] + 0.1, centralityBins[i+1] - 0.1);
2647         TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
2648         if (verbose)
2649         {
2650           new TCanvas;
2651           proj->DrawCopy();
2652         }
2653         if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
2654           target[i] = proj->GetFunction("pol0")->GetParameter(0);
2655         else
2656           target[i] = 0;
2657         Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
2658       }
2659     }
2660
2661     // extend for full pT range
2662     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2663       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2664         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2665         {
2666           for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
2667           {
2668
2669             Int_t bins[4];
2670             bins[0] = x;
2671             bins[1] = y;
2672             bins[2] = z;
2673             bins[3] = z2;
2674
2675             Int_t z2Bin = 0;
2676             while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
2677               z2Bin++;
2678
2679             //Printf("%d %d", z2, z2Bin);
2680
2681             fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2682             fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
2683             if (trackingCont[0] > 0)
2684               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
2685             else
2686               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
2687           }
2688        }
2689
2690      delete[] trackingEff;
2691      delete[] trackingCont;
2692    }
2693
2694    SetCentralityRange(0, 0);
2695 }
2696
2697 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
2698 {
2699   // corrects the dphi distribution with an extra factor close to dphi ~ 0
2700
2701   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
2702
2703   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
2704
2705   // optimized implementation
2706   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2707   {
2708     Int_t bins[6];
2709     Double_t value = grid->GetBinContent(binIdx, bins);
2710     Double_t error = grid->GetBinError(binIdx);
2711
2712     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
2713     if (TMath::Abs(binCenter) < 0.2)
2714     {
2715       value *= 0.985;
2716       error *= 0.985;
2717     }
2718     else if (TMath::Abs(binCenter) < 0.3)
2719     {
2720       value *= 0.9925;
2721       error *= 0.9925;
2722     }
2723
2724     grid->SetBinContent(bins, value);
2725     grid->SetBinError(bins, error);
2726   }
2727 }
2728
2729 void AliUEHist::Scale(Double_t factor)
2730 {
2731   // scales all contained histograms by the given factor
2732
2733   for (Int_t i=0; i<4; i++)
2734     if (fTrackHist[i])
2735       fTrackHist[i]->Scale(factor);
2736
2737   fEventHist->Scale(factor);
2738   fTrackHistEfficiency->Scale(factor);
2739 }
2740
2741 void AliUEHist::Reset()
2742 {
2743   // resets all contained histograms
2744
2745   for (Int_t i=0; i<4; i++)
2746     if (fTrackHist[i])
2747       for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
2748         fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2749
2750   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2751     fEventHist->GetGrid(step)->GetGrid()->Reset();
2752
2753   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
2754     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
2755 }
2756
2757 THnBase* AliUEHist::ChangeToThn(THnBase* sparse)
2758 {
2759   // change the object to THn for faster processing
2760
2761         // convert to THn (SEGV's for some strange reason...)
2762         // x = THn::CreateHn("a", "a", sparse);
2763
2764   // own implementation
2765   Int_t nBins[10];
2766   for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2767     nBins[i] = sparse->GetAxis(i)->GetNbins();
2768   THn* tmpTHn = new THnF(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse->GetNdimensions(), nBins, 0, 0);
2769   for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2770   {
2771     tmpTHn->SetBinEdges(i, sparse->GetAxis(i)->GetXbins()->GetArray());
2772     tmpTHn->GetAxis(i)->SetTitle(sparse->GetAxis(i)->GetTitle());
2773   }
2774   tmpTHn->RebinnedAdd(sparse);
2775
2776   return tmpTHn;
2777 }
2778
2779 void AliUEHist::CondenseBin(THnSparse* grid, THnSparse* target, Int_t axis, Float_t targetValue, Float_t from, Float_t to)
2780 {
2781   //
2782   // loops through the histogram and moves all entries to a single point <targetValue> on the axis <axis>
2783   // if <from> and <to> are set, then moving only occurs if value on <axis> is betweem <from> and <to>
2784   //
2785
2786   if (grid->GetNdimensions() > 6)
2787     AliFatal("Too many dimensions in THnSparse");
2788
2789   Int_t targetBin = grid->GetAxis(axis)->FindBin(targetValue);
2790   AliInfo(Form("Target bin on axis %d with value %f is %d", axis, targetValue, targetBin));
2791
2792   Int_t fromBin = 1;
2793   Int_t toBin = grid->GetAxis(axis)->GetNbins();
2794   if (to > from)
2795   {
2796     fromBin = grid->GetAxis(axis)->FindBin(from);
2797     toBin = grid->GetAxis(axis)->FindBin(to);
2798     AliInfo(Form("Only condensing between bin %d and %d", fromBin, toBin));
2799   }
2800
2801   Int_t bins[6];
2802   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2803   {
2804     Double_t value = grid->GetBinContent(binIdx, bins);
2805     Double_t error = grid->GetBinError(binIdx);
2806
2807     if (bins[axis] >= fromBin && bins[axis] <= toBin)
2808       bins[axis] = targetBin;
2809
2810     value += target->GetBinContent(bins);
2811     error = TMath::Sqrt(error * error + target->GetBinError(bins) * target->GetBinError(bins));
2812
2813     target->SetBinContent(bins, value);
2814     target->SetBinError(bins, error);
2815   }
2816 }
2817
2818 void AliUEHist::CondenseBin(CFStep step, Int_t trackAxis, Int_t eventAxis, Float_t targetValue, Float_t from, Float_t to, CFStep tmpStep)
2819 {
2820   // loops through the histogram at <step> and moves all entries to a single point <targetValue> on the axes
2821   // <trackAxis> and <eventAxis>. <tmpStep> is used to temporary store the data
2822   // This is useful e.g. to move bin content around for MC productions where the centrality selection did
2823   // not yield the desired result
2824
2825   // reset tmpStep
2826   fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2827   for (UInt_t i=0; i<fkRegions; i++)
2828     if (fTrackHist[i])
2829       fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2830
2831   // copy to tmpStep
2832   CorrectTracks(step, tmpStep, 0, -1);
2833   CorrectEvents(step, tmpStep, 0, -1);
2834
2835   // reset step
2836   fEventHist->GetGrid(step)->GetGrid()->Reset();
2837   for (UInt_t i=0; i<fkRegions; i++)
2838     if (fTrackHist[i])
2839       fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2840
2841   // rewriting
2842   for (UInt_t i=0; i<fkRegions; i++)
2843   {
2844     if (!fTrackHist[i])
2845       continue;
2846
2847     THnSparse* grid = fTrackHist[i]->GetGrid(tmpStep)->GetGrid();
2848     THnSparse* target = fTrackHist[i]->GetGrid(step)->GetGrid();
2849
2850     CondenseBin(grid, target, trackAxis, targetValue, from, to);
2851   }
2852   CondenseBin(fEventHist->GetGrid(tmpStep)->GetGrid(), fEventHist->GetGrid(step)->GetGrid(), eventAxis, targetValue, from, to);
2853
2854   // reset tmpStep
2855   fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2856   for (UInt_t i=0; i<fkRegions; i++)
2857     if (fTrackHist[i])
2858       fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2859 }