PWG4/JetTasks/AliUEHist.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliUEHist.cxx 20164 2007-08-14 15:31:50Z morsch $ */
  17
  18 //
  19 //
  20 // encapsulate histogram and corrections for one underlying event histogram
  21 //
  22 //
  23 // Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero
  24
  25 #include "AliUEHist.h"
  26 #include "AliCFContainer.h"
  27 #include "THnSparse.h"
  28 #include "TMath.h"
  29 #include "TList.h"
  30 #include "TCollection.h"
  31 #include "TH1D.h"
  32 #include "TH2D.h"
  33 #include "TH3D.h"
  34 #include "AliLog.h"
  35 #include "TCanvas.h"
  36 #include "TF1.h"
  37 #include "AliTHn.h"
  38
  39 ClassImp(AliUEHist)
  40
  41 const Int_t AliUEHist::fgkCFSteps = 10;
  42
  43 AliUEHist::AliUEHist(const char* reqHist) :
  44   TObject(),
  45   fkRegions(4),
  46   fEventHist(0),
  47   fTrackHistEfficiency(0),
  48   fEtaMin(0),
  49   fEtaMax(0),
  50   fPtMin(0),
  51   fPtMax(0),
  52   fCentralityMin(0),
  53   fCentralityMax(0),
  54   fZVtxMin(0),
  55   fZVtxMax(0),
  56   fContaminationEnhancement(0),
  57   fCombineMinMax(0),
  58   fCache(0),
  59   fHistogramType(reqHist)
  60 {
  61   // Constructor
  62
  63   for (UInt_t i=0; i<fkRegions; i++)
  64     fTrackHist[i] = 0;
  65
  66   if (strlen(reqHist) == 0)
  67     return;
  68
  69   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  70   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  71
  72   const char* title = "";
  73
  74   // track level
  75   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  76   Int_t iTrackBin[6];
  77   Double_t* trackBins[6];
  78   const char* trackAxisTitle[6];
  79
  80   // eta
  81   iTrackBin[0] = 20;
  82   Double_t etaBins[20+1];
  83   for (Int_t i=0; i<=iTrackBin[0]; i++)
  84     etaBins[i] = -1.0 + 0.1 * i;
  85   trackBins[0] = etaBins;
  86   trackAxisTitle[0] = "#eta";
  87
  88   // delta eta
  89   const Int_t kNDeltaEtaBins = 40;
  90   Double_t deltaEtaBins[kNDeltaEtaBins+1];
  91   for (Int_t i=0; i<=kNDeltaEtaBins; i++)
  92     deltaEtaBins[i] = -2.0 + 0.1 * i;
  93
  94   // pT
  95   iTrackBin[1] = 22;
  96   Double_t pTBins[] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0, 20.0};
  97   trackBins[1] = pTBins;
  98   trackAxisTitle[1] = "p_{T} (GeV/c)";
  99
 100   // pT,lead binning 1
 101   const Int_t kNLeadingpTBins = 100;
 102   Double_t leadingpTBins[kNLeadingpTBins+1];
 103   for (Int_t i=0; i<=kNLeadingpTBins; i++)
 104     leadingpTBins[i] = 0.5 * i;
 105
 106   // pT,lead binning 2
 107   const Int_t kNLeadingpTBins2 = 9;
 108   Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0 };
 109
 110   // phi,lead; this binning starts at -pi/2 and is modulo 3
 111   const Int_t kNLeadingPhiBins = 72;
 112   Double_t leadingPhiBins[kNLeadingPhiBins+1];
 113   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
 114     leadingPhiBins[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiBins * i * TMath::TwoPi();
 115
 116   // multiplicity
 117   const Int_t kNMultiplicityBins = 15;
 118   Double_t multiplicityBins[kNMultiplicityBins+1];
 119   for (Int_t i=0; i<=kNMultiplicityBins; i++)
 120     multiplicityBins[i] = -0.5 + i;
 121   multiplicityBins[kNMultiplicityBins] = 200;
 122
 123   trackBins[3] = multiplicityBins;
 124   iTrackBin[3] = kNMultiplicityBins;
 125   trackAxisTitle[3] = "multiplicity";
 126
 127   // centrality (in %)
 128   const Int_t kNCentralityBins = 5 + 1 + 9;
 129   Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
 130
 131   // particle species
 132   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
 133   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
 134
 135   // vtx-z axis
 136   const Int_t kNVertexBins = 7;
 137   Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
 138
 139   Bool_t useVtxAxis = kFALSE;
 140
 141   // selection depending on requested histogram
 142   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
 143   if (strcmp(reqHist, "NumberDensitypT") == 0)
 144   {
 145     axis = 0;
 146     title = "d^{2}N_{ch}/d#varphid#eta";
 147   }
 148   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 149   {
 150     axis = 1;
 151     title = "d^{2}N_{ch}/d#varphid#eta";
 152   }
 153   else if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0 || strcmp(reqHist, "NumberDensityPhiCentralityVtx") == 0)
 154   {
 155     if (strcmp(reqHist, "NumberDensityPhiCentralityVtx") == 0)
 156     {
 157       reqHist = "NumberDensityPhiCentrality";
 158       fHistogramType = reqHist;
 159       useVtxAxis = kTRUE;
 160     }
 161     axis = 2;
 162     title = "d^{2}N_{ch}/d#varphid#eta";
 163   }
 164   else if (strcmp(reqHist, "SumpT") == 0)
 165   {
 166     axis = 0;
 167     title = "d^{2}#Sigma p_{T}/d#varphid#eta";
 168   }
 169   else
 170     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 171
 172   UInt_t initRegions = fkRegions;
 173
 174   if (axis == 0)
 175   {
 176     trackBins[2] = leadingpTBins;
 177     iTrackBin[2] = kNLeadingpTBins;
 178     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 179
 180   }
 181   else if (axis == 1)
 182   {
 183     nTrackVars = 5;
 184     initRegions = 1;
 185
 186     iTrackBin[2] = kNLeadingpTBins2;
 187     trackBins[2] = leadingpTBins2;
 188     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 189
 190     iTrackBin[4] = kNLeadingPhiBins;
 191     trackBins[4] = leadingPhiBins;
 192     trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
 193   }
 194   else if (axis == 2)
 195   {
 196     nTrackVars = 5;
 197     initRegions = 1;
 198
 199     iTrackBin[0] = kNDeltaEtaBins;
 200     trackBins[0] = deltaEtaBins;
 201     trackAxisTitle[0] = "#Delta#eta";
 202
 203     iTrackBin[2] = kNLeadingpTBins2;
 204     trackBins[2] = leadingpTBins2;
 205     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 206
 207     trackBins[3] = centralityBins;
 208     iTrackBin[3] = kNCentralityBins;
 209     trackAxisTitle[3] = "centrality";
 210
 211     iTrackBin[4] = kNLeadingPhiBins;
 212     trackBins[4] = leadingPhiBins;
 213     trackAxisTitle[4] = "#Delta#varphi (rad.)";
 214
 215     if (useVtxAxis)
 216     {
 217       nTrackVars = 6;
 218       iTrackBin[5] = kNVertexBins;
 219       trackBins[5] = vertexBins;
 220       trackAxisTitle[5] = "z-vtx (cm)";
 221     }
 222   }
 223
 224   for (UInt_t i=0; i<initRegions; i++)
 225   {
 226     if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
 227       fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 228     else
 229       fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 230
 231     for (Int_t j=0; j<nTrackVars; j++)
 232     {
 233       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 234       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 235     }
 236
 237     SetStepNames(fTrackHist[i]);
 238   }
 239
 240   // event level
 241   Int_t nEventVars = 2;
 242   Int_t iEventBin[3];
 243
 244   // track 3rd and 4th axis --> event 1st and 2nd axis
 245   iEventBin[0] = iTrackBin[2];
 246   iEventBin[1] = iTrackBin[3];
 247
 248   // plus track 5th axis (in certain cases)
 249   if (axis == 2 && useVtxAxis)
 250   {
 251     nEventVars = 3;
 252     iEventBin[2] = iTrackBin[5];
 253   }
 254
 255   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
 256
 257   fEventHist->SetBinLimits(0, trackBins[2]);
 258   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 259
 260   fEventHist->SetBinLimits(1, trackBins[3]);
 261   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 262
 263   if (axis == 2 && useVtxAxis)
 264   {
 265     fEventHist->SetBinLimits(2, trackBins[5]);
 266     fEventHist->SetVarTitle(2, trackAxisTitle[5]);
 267   }
 268
 269   SetStepNames(fEventHist);
 270
 271   iTrackBin[2] = kNSpeciesBins;
 272
 273   fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 3, 4, iTrackBin);
 274   fTrackHistEfficiency->SetBinLimits(0, trackBins[0]);
 275   fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]);
 276   fTrackHistEfficiency->SetBinLimits(1, trackBins[1]);
 277   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
 278   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
 279   fTrackHistEfficiency->SetVarTitle(2, "particle species");
 280   fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
 281   fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
 282 }
 283
 284 //_____________________________________________________________________________
 285 AliUEHist::AliUEHist(const AliUEHist &c) :
 286   TObject(),
 287   fkRegions(4),
 288   fEventHist(0),
 289   fTrackHistEfficiency(0),
 290   fEtaMin(0),
 291   fEtaMax(0),
 292   fPtMin(0),
 293   fPtMax(0),
 294   fCentralityMin(0),
 295   fCentralityMax(0),
 296   fZVtxMin(0),
 297   fZVtxMax(0),
 298   fContaminationEnhancement(0),
 299   fCombineMinMax(0),
 300   fCache(0),
 301   fHistogramType()
 302 {
 303   //
 304   // AliUEHist copy constructor
 305   //
 306
 307   ((AliUEHist &) c).Copy(*this);
 308 }
 309
 310 //____________________________________________________________________
 311 void AliUEHist::SetStepNames(AliCFContainer* container)
 312 {
 313   // sets the names of the correction steps
 314
 315   for (Int_t i=0; i<fgkCFSteps; i++)
 316     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 317 }
 318
 319 //____________________________________________________________________
 320 AliUEHist::~AliUEHist()
 321 {
 322   // Destructor
 323
 324   for (UInt_t i=0; i<fkRegions; i++)
 325   {
 326     if (fTrackHist[i])
 327     {
 328       delete fTrackHist[i];
 329       fTrackHist[i] = 0;
 330     }
 331   }
 332
 333   if (fEventHist)
 334   {
 335     delete fEventHist;
 336     fEventHist = 0;
 337   }
 338
 339   if (fTrackHistEfficiency)
 340   {
 341     delete fTrackHistEfficiency;
 342     fTrackHistEfficiency = 0;
 343   }
 344
 345   if (fCache)
 346   {
 347     delete fCache;
 348     fCache = 0;
 349   }
 350 }
 351
 352 //____________________________________________________________________
 353 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 354 {
 355   // assigment operator
 356
 357   if (this != &c)
 358     ((AliUEHist &) c).Copy(*this);
 359
 360   return *this;
 361 }
 362
 363 //____________________________________________________________________
 364 void AliUEHist::Copy(TObject& c) const
 365 {
 366   // copy function
 367
 368   AliUEHist& target = (AliUEHist &) c;
 369
 370   for (UInt_t i=0; i<fkRegions; i++)
 371     if (fTrackHist[i])
 372       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 373
 374   if (fEventHist)
 375     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 376
 377   if (fTrackHistEfficiency)
 378     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
 379
 380   target.fEtaMin = fEtaMin;
 381   target.fEtaMax = fEtaMax;
 382   target.fPtMin = fPtMin;
 383   target.fPtMax = fPtMax;
 384   target.fCentralityMin = fCentralityMin;
 385   target.fCentralityMax = fCentralityMax;
 386   target.fZVtxMin = fZVtxMin;
 387   target.fZVtxMax = fZVtxMax;
 388
 389   if (fContaminationEnhancement)
 390     target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
 391
 392   target.fCombineMinMax = fCombineMinMax;
 393   target.fHistogramType = fHistogramType;
 394 }
 395
 396 //____________________________________________________________________
 397 Long64_t AliUEHist::Merge(TCollection* list)
 398 {
 399   // Merge a list of AliUEHist objects with this (needed for
 400   // PROOF).
 401   // Returns the number of merged objects (including this).
 402
 403   if (!list)
 404     return 0;
 405
 406   if (list->IsEmpty())
 407     return 1;
 408
 409   TIterator* iter = list->MakeIterator();
 410   TObject* obj;
 411
 412   // collections of objects
 413   const UInt_t kMaxLists = fkRegions+2;
 414   TList** lists = new TList*[kMaxLists];
 415
 416   for (UInt_t i=0; i<kMaxLists; i++)
 417     lists[i] = new TList;
 418
 419   Int_t count = 0;
 420   while ((obj = iter->Next())) {
 421
 422     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 423     if (entry == 0)
 424       continue;
 425
 426     for (UInt_t i=0; i<fkRegions; i++)
 427       if (entry->fTrackHist[i])
 428         lists[i]->Add(entry->fTrackHist[i]);
 429
 430     lists[fkRegions]->Add(entry->fEventHist);
 431     lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
 432
 433     count++;
 434   }
 435   for (UInt_t i=0; i<fkRegions; i++)
 436     if (fTrackHist[i])
 437       fTrackHist[i]->Merge(lists[i]);
 438
 439   fEventHist->Merge(lists[fkRegions]);
 440   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
 441
 442   for (UInt_t i=0; i<kMaxLists; i++)
 443     delete lists[i];
 444
 445   delete[] lists;
 446
 447   return count+1;
 448 }
 449
 450 //____________________________________________________________________
 451 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 452 {
 453   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 454
 455   if (fEtaMax > fEtaMin)
 456     grid->SetRangeUser(0, fEtaMin, fEtaMax);
 457   if (fPtMax > fPtMin)
 458     grid->SetRangeUser(1, fPtMin, fPtMax);
 459 }
 460
 461 //____________________________________________________________________
 462 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 463 {
 464   // resets all bin limits
 465
 466   for (Int_t i=0; i<grid->GetNVar(); i++)
 467     if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange))
 468       grid->SetRangeUser(i, 0, -1);
 469 }
 470
 471 //____________________________________________________________________
 472 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 473 {
 474   // prints the number of empty bins in the track end event histograms in the given step
 475
 476   Int_t binBegin[4];
 477   Int_t binEnd[4];
 478
 479   for (Int_t i=0; i<4; i++)
 480   {
 481     binBegin[i] = 1;
 482     binEnd[i]   = fTrackHist[0]->GetNBins(i);
 483   }
 484
 485   if (fEtaMax > fEtaMin)
 486   {
 487     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
 488     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
 489   }
 490
 491   if (fPtMax > fPtMin)
 492   {
 493     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
 494     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
 495   }
 496
 497   if (ptLeadMax > ptLeadMin)
 498   {
 499     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
 500     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
 501   }
 502
 503   // start from multiplicity 1
 504   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
 505
 506   for (UInt_t region=0; region<fkRegions; region++)
 507   {
 508     Int_t total = 0;
 509     Int_t count = 0;
 510     Int_t vars[4];
 511
 512     for (Int_t i=0; i<4; i++)
 513       vars[i] = binBegin[i];
 514
 515     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
 516     while (1)
 517     {
 518       if (grid->GetElement(vars) == 0)
 519       {
 520         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
 521           grid->GetBinCenter(0, vars[0]),
 522           grid->GetBinCenter(1, vars[1]),
 523           grid->GetBinCenter(2, vars[2]),
 524           grid->GetBinCenter(3, vars[3])
 525         );
 526         count++;
 527       }
 528
 529       vars[3]++;
 530       for (Int_t i=3; i>0; i--)
 531       {
 532         if (vars[i] == binEnd[i]+1)
 533         {
 534           vars[i] = binBegin[i];
 535           vars[i-1]++;
 536         }
 537       }
 538
 539       if (vars[0] == binEnd[0]+1)
 540         break;
 541       total++;
 542     }
 543
 544     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
 545   }
 546 }
 547
 548 //____________________________________________________________________
 549 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)
 550 {
 551   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 552   //
 553   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 554   // Histogram has to be deleted by the caller of the function
 555   //
 556   // twoD: 0: 1D histogram as function of phi
 557   //       1: 2D histogram, deltaphi, deltaeta
 558   //       10: 1D histogram, within |deltaeta| < 0.8
 559   //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
 560   //
 561   // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
 562   //
 563   // normEvents: if non-0 the number of events/trigger particles for the normalization is filled
 564
 565   // unzoom all axes
 566   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 567   ResetBinLimits(fEventHist->GetGrid(step));
 568
 569   SetBinLimits(fTrackHist[region]->GetGrid(step));
 570
 571   // z vtx
 572   if (fZVtxMax > fZVtxMin)
 573   {
 574     Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
 575     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
 576     fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
 577   }
 578
 579   TH1D* tracks = 0;
 580
 581   if (ptLeadMin < 0)
 582   {
 583     tracks = fTrackHist[region]->ShowProjection(2, step);
 584     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 585     if (fCombineMinMax && region == kMin)
 586     {
 587       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 588       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 589
 590       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 591       tracks->Add(tracks2);
 592
 593       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 594     }
 595
 596     // normalize to get a density (deta dphi)
 597     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 598     Float_t phiRegion = TMath::TwoPi() / 3;
 599     if (!fCombineMinMax && region == kMin)
 600       phiRegion /= 2;
 601     Float_t normalization = phiRegion;
 602     if (etaNorm)
 603       normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 604     //Printf("Normalization: %f", normalization);
 605     tracks->Scale(1.0 / normalization);
 606
 607     TH1D* events = fEventHist->ShowProjection(0, step);
 608     tracks->Divide(events);
 609   }
 610   else
 611   {
 612     if (multBinEnd >= multBinBegin)
 613     {
 614       Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 615       fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 616       fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 617     }
 618
 619     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 620     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 621
 622     Printf("Using leading pT range %d --> %d", firstBin, lastBin);
 623
 624     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
 625
 626     if (twoD == 0)
 627       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
 628     else
 629       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
 630
 631     Printf("Calculated histogram --> %f tracks", tracks->Integral());
 632     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 633
 634     if (twoD == 10 || twoD == 11)
 635     {
 636       Float_t etaLimit = 1.0;
 637       if (twoD == 10)
 638       {
 639         tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
 640
 641         // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
 642         tracks->Scale(1. / 0.75);
 643       }
 644       else if (twoD == 11)
 645       {
 646         TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
 647         TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
 648
 649         tracksTmp1->Add(tracksTmp2);
 650
 651         delete tracks;
 652         tracks = tracksTmp1;
 653         delete tracksTmp2;
 654
 655         tracks->Scale(1. / 0.25);
 656       }
 657     }
 658
 659     // normalize to get a density (deta dphi)
 660     // TODO normalization may be off for 2d histogram
 661     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
 662
 663     if (etaNorm)
 664     {
 665       TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 666       if (strcmp(axis->GetTitle(), "#eta") == 0)
 667       {
 668         Printf("Normalizing using eta-axis range");
 669         normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
 670       }
 671       else
 672       {
 673         Printf("Normalizing assuming accepted range of |eta| < 0.8");
 674         normalization *= 0.8 * 2;
 675       }
 676     }
 677
 678     //Printf("Normalization: %f", normalization);
 679     tracks->Scale(1.0 / normalization);
 680
 681     // 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!
 682     TH1D* events = fEventHist->ShowProjection(0, step);
 683     Int_t nEvents = (Int_t) events->Integral(firstBin, lastBin);
 684     Printf("Calculated histogram --> %d events", nEvents);
 685     if (normEvents)
 686       *normEvents = nEvents;
 687
 688     if (nEvents > 0)
 689       tracks->Scale(1.0 / nEvents);
 690
 691     delete events;
 692   }
 693
 694   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 695   ResetBinLimits(fEventHist->GetGrid(step));
 696
 697   return tracks;
 698 }
 699
 700 //____________________________________________________________________
 701 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)
 702 {
 703   // Calls GetUEHist(...) for *each* multiplicity bin and performs a sum of ratios:
 704   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
 705   // where N is the total number of events/trigger particles and the subscript is the multiplicity bin
 706   //
 707   // Can only be used for the 2D histogram at present
 708   //
 709   // Parameters:
 710   //   mixed: AliUEHist containing mixed event corresponding to this object
 711   //   <other parameters> : check documentation of AliUEHist::GetUEHist
 712
 713   Int_t multIter = multBinBegin;
 714
 715   TH2* totalTracks = 0;
 716   Int_t totalEvents = 0;
 717
 718   Int_t vertexBin = 1;
 719   TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
 720   if (useVertexBins && !vertexAxis)
 721   {
 722     Printf("Vertex axis requested but not available");
 723     return 0;
 724   }
 725
 726   // vertex bin loop
 727   while (1)
 728   {
 729     if (useVertexBins)
 730     {
 731       SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
 732       mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
 733       vertexBin++;
 734     }
 735
 736     // multiplicity loop
 737     while (1)
 738     {
 739       Int_t multBinBeginLocal = multBinBegin;
 740       Int_t multBinEndLocal = multBinEnd;
 741
 742       if (multBinEnd >= multBinBegin)
 743       {
 744         multBinBeginLocal = multIter;
 745         multBinEndLocal = multIter;
 746         multIter++;
 747       }
 748
 749       Int_t nEvents = 0;
 750       TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
 751       // undo normalization
 752       tracks->Scale(nEvents);
 753       totalEvents += nEvents;
 754
 755       TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
 756
 757       // asssume flat in dphi, gain in statistics
 758   //     TH1* histMixedproj = mixedTwoD->ProjectionY();
 759   //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
 760   //
 761   //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
 762   //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
 763   //    mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
 764
 765   //       delete histMixedproj;
 766
 767       // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
 768 /*      Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
 769       mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
 770       tracks->Scale(mixedNorm);*/
 771
 772       tracks->Scale(mixedTwoD->Integral() / tracks->Integral());
 773
 774       tracks->Divide(mixedTwoD);
 775
 776       delete mixedTwoD;
 777
 778       if (!totalTracks)
 779         totalTracks = tracks;
 780       else
 781       {
 782         totalTracks->Add(tracks);
 783         delete tracks;
 784       }
 785
 786       if (multIter > multBinEnd)
 787         break;
 788     }
 789
 790     if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
 791       break;
 792   }
 793
 794   if (useVertexBins)
 795     totalEvents = vertexAxis->GetNbins();
 796   Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
 797   if (totalEvents > 0)
 798     totalTracks->Scale(1.0 / totalEvents);
 799
 800   return totalTracks;
 801 }
 802
 803 //____________________________________________________________________
 804 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)
 805 {
 806   // returns a histogram projected to pT,assoc with the provided cuts
 807
 808    // unzoom all axes
 809   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 810   ResetBinLimits(fEventHist->GetGrid(step));
 811
 812   TH1D* tracks = 0;
 813
 814     // 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
 815   // therefore the number density must be calculated here per leading pT bin and then summed
 816
 817   if (multBinEnd > multBinBegin)
 818     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
 819   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
 820   fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
 821
 822   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
 823   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
 824
 825   // get real phi cuts due to binning
 826   Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
 827   Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
 828   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);
 829
 830   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
 831   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
 832
 833   TH1D* events = fEventHist->ShowProjection(0, step);
 834
 835   for (Int_t bin=firstBin; bin<=lastBin; bin++)
 836   {
 837     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
 838
 839     // project to pT,assoc
 840     TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
 841
 842     Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
 843     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 844
 845     // normalize to get a density (deta dphi)
 846     Float_t normalization = 1;
 847     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 848     if (strcmp(axis->GetTitle(), "#eta") == 0)
 849     {
 850       Printf("Normalizing using eta-axis range");
 851       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
 852     }
 853     else
 854     {
 855       Printf("Normalizating assuming accepted range of |eta| < 0.8");
 856       normalization *= 0.8 * 2;
 857     }
 858
 859     // dphi
 860     if (!skipPhiNormalization)
 861       normalization *= phiMaxReal - phiMinReal;
 862
 863     //Printf("Normalization: %f", normalization);
 864     tracksTmp->Scale(1.0 / normalization);
 865
 866     // and now dpT (bins have different width)
 867     for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
 868     {
 869       tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
 870       tracksTmp->SetBinError  (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
 871     }
 872
 873     Int_t nEvents = (Int_t) events->GetBinContent(bin);
 874     if (nEvents > 0)
 875       tracksTmp->Scale(1.0 / nEvents);
 876     Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
 877
 878     if (!tracks)
 879       tracks = tracksTmp;
 880     else
 881     {
 882       tracks->Add(tracksTmp);
 883       delete tracksTmp;
 884     }
 885   }
 886
 887   delete events;
 888
 889   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 890   ResetBinLimits(fEventHist->GetGrid(step));
 891
 892   return tracks;
 893 }
 894
 895 void AliUEHist::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
 896 {
 897   // multiplies <grid> with <histogram> and put the result in <target>
 898   // <grid> has usually more dimensions than histogram. The axis which are used to choose the value
 899   // from <histogram> are given in <var1> and <var2>
 900   //
 901   // if <histogram> is 0, just copies content from step1 to step2
 902
 903   // clear target histogram
 904   target->Reset();
 905
 906   if (histogram != 0)
 907   {
 908     if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
 909       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
 910
 911     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
 912       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
 913   }
 914
 915   if (grid->GetNdimensions() > 6)
 916     AliFatal("Too many dimensions in THnSparse");
 917
 918   Int_t bins[6];
 919
 920   // optimized implementation
 921   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 922   {
 923     Double_t value = grid->GetBinContent(binIdx, bins);
 924     Double_t error = grid->GetBinError(binIdx);
 925
 926     if (histogram != 0)
 927     {
 928       if (var2 < 0)
 929       {
 930         value *= histogram->GetBinContent(bins[var1]);
 931         error *= histogram->GetBinContent(bins[var1]);
 932       }
 933       else
 934       {
 935         value *= histogram->GetBinContent(bins[var1], bins[var2]);
 936         error *= histogram->GetBinContent(bins[var1], bins[var2]);
 937       }
 938     }
 939
 940     target->SetBinContent(bins, value);
 941     target->SetBinError(bins, error);
 942   }
 943 }
 944
 945 //____________________________________________________________________
 946 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 947 {
 948   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
 949   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
 950   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
 951   //
 952   // if trackCorrection is 0, just copies content from step1 to step2
 953
 954   for (UInt_t region=0; region<fkRegions; region++)
 955     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
 956 }
 957
 958 //____________________________________________________________________
 959 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
 960 {
 961   //
 962   // see documentation of CorrectTracks above
 963   //
 964
 965   if (!fTrackHist[region])
 966     return;
 967
 968   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
 969   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
 970
 971   MultiplyHistograms(grid, target, trackCorrection, var1, var2);
 972
 973   Printf("AliUEHist::CorrectTracks: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0);
 974 }
 975
 976 //____________________________________________________________________
 977 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
 978 {
 979   // corrects from step1 to step2 by multiplying the events with eventCorrection
 980   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
 981   //
 982   // if eventCorrection is 0, just copies content from step1 to step2
 983
 984   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
 985   AliCFGridSparse* target = fEventHist->GetGrid(step2);
 986
 987   MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);
 988
 989   Printf("AliUEHist::CorrectEvents: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0);
 990 }
 991
 992 //____________________________________________________________________
 993 void AliUEHist::Correct(AliUEHist* corrections)
 994 {
 995   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
 996   //
 997   // in this object the data is expected in the step kCFStepReconstructed
 998
 999   if (fHistogramType.Length() == 0)
1000   {
1001     Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
1002     if (fTrackHist[kToward]->GetNVar() < 5)
1003     {
1004       if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#varphid#eta") == 0)
1005         fHistogramType = "NumberDensitypT";
1006       else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#eta") == 0)
1007         fHistogramType = "SumpT";
1008     }
1009     else if (fTrackHist[kToward]->GetNVar() == 5)
1010     {
1011       if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
1012         fHistogramType = "NumberDensityPhi";
1013       else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
1014         fHistogramType = "NumberDensityPhiCentrality";
1015     }
1016
1017     if (fHistogramType.Length() == 0)
1018       AliFatal("Cannot figure out histogram type. Try setting it manually...");
1019   }
1020
1021   Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
1022
1023   if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
1024   {
1025     // ---- track level
1026
1027     // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
1028     // extracted as function of leading pT
1029     Bool_t biasFromMC = kFALSE;
1030     const char* projAxis = "z";
1031     Int_t secondBin = -1;
1032
1033     if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
1034     {
1035       projAxis = "yz";
1036       secondBin = 4;
1037     }
1038
1039     for (UInt_t region = 0; region < fkRegions; region++)
1040     {
1041       if (!fTrackHist[region])
1042         continue;
1043
1044       TH1* leadingBiasTracks =  0;
1045       if (biasFromMC)
1046       {
1047         leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
1048         Printf("WARNING: Using MC bias correction");
1049       }
1050       else
1051         leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
1052
1053       CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
1054       if (region == kMin && fCombineMinMax)
1055       {
1056         CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
1057         delete leadingBiasTracks;
1058         break;
1059       }
1060       delete leadingBiasTracks;
1061     }
1062
1063     TH1* leadingBiasEvents = 0;
1064     if (biasFromMC)
1065       leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
1066     else
1067       leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
1068
1069     //new TCanvas; leadingBiasEvents->DrawCopy();
1070     CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
1071     delete leadingBiasEvents;
1072
1073     // --- contamination correction ---
1074
1075     TH2D* contamination = corrections->GetTrackingContamination();
1076     if (corrections->fContaminationEnhancement)
1077     {
1078       Printf("Applying contamination enhancement");
1079
1080       for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
1081         for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
1082           contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
1083     }
1084     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
1085     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1086     delete contamination;
1087
1088     // --- efficiency correction ---
1089     // correct single-particle efficiency for associated particles
1090     // in addition correct for efficiency on trigger particles (tracks AND events)
1091
1092     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
1093     // use kCFStepVertex as a temporary step (will be overwritten later anyway)
1094     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1095     delete efficiencyCorrection;
1096
1097     // correct pT,T
1098     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
1099     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
1100     CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
1101     delete efficiencyCorrection;
1102
1103     // copy
1104     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
1105     CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
1106
1107     // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
1108     // practically independent of low pT cut
1109     TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
1110
1111     // convert stage from true multiplicity to observed multiplicity by simple conversion factor
1112     TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
1113     vertexCorrectionObs->Reset();
1114
1115     TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
1116     // some defaults
1117     func->SetParameters(0.1, 1, -0.7);
1118     vertexCorrection->Fit(func, "0I", "", 0, 3);
1119     for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
1120     {
1121       Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
1122       if (xPos < 3)
1123         vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
1124       else
1125         vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
1126     }
1127
1128     #if 0
1129     new TCanvas;
1130     vertexCorrection->DrawCopy();
1131     vertexCorrectionObs->SetLineColor(2);
1132     vertexCorrectionObs->DrawCopy("same");
1133     func->SetRange(0, 4);
1134     func->DrawClone("same");
1135     #endif
1136
1137     CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
1138     CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
1139     delete vertexCorrectionObs;
1140     delete vertexCorrection;
1141     delete func;
1142
1143     // copy
1144     CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
1145     CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
1146   }
1147   else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
1148   {
1149     // copy
1150     CorrectTracks(kCFStepReconstructed, kCFStepTracked, 0, -1);
1151     CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, -1);
1152
1153     // Dont use eta in the following, because it is a Delta-eta axis
1154
1155     // contamination correction
1156     // correct single-particle contamination for associated particles
1157
1158     TH1* contamination = corrections->GetTrackingContamination(1);
1159
1160     if (1)
1161     {
1162       Printf("Applying contamination enhancement");
1163
1164       for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1165       {
1166         printf("%f", contamination->GetBinContent(bin));
1167         if (contamination->GetBinContent(bin) > 0)
1168           contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1169         printf(" --> %f\n", contamination->GetBinContent(bin));
1170       }
1171     }
1172
1173     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1174     delete contamination;
1175
1176     // correct for additional contamination due to trigger particle around phi ~ 0
1177     TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1178     if (1)
1179     {
1180       Printf("Applying contamination enhancement");
1181
1182       for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1183         for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1184         {
1185           printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1186           if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1187             correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1188           printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1189         }
1190     }
1191
1192     new TCanvas; correlatedContamination->DrawCopy("COLZ");
1193     CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1194 //     Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1195
1196     delete correlatedContamination;
1197
1198     // TODO correct for contamination of trigger particles (for tracks AND events)
1199     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1200
1201     // --- efficiency correction ---
1202     // correct single-particle efficiency for associated particles
1203     // in addition correct for efficiency on trigger particles (tracks AND events)
1204
1205     // in bins of pT and centrality
1206     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1207     // use kCFStepAnaTopology as a temporary step
1208     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1209     delete efficiencyCorrection;
1210
1211     // correct pT,T in bins of pT and centrality
1212     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1213     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1214     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1215     delete efficiencyCorrection;
1216
1217     // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1218     // copy
1219     CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1220     CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1221   }
1222   else
1223     AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1224 }
1225
1226 //____________________________________________________________________
1227 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1228 {
1229   // creates a track-level efficiency by dividing step2 by step1
1230   // projected to axis1 and axis2 (optional if >= 0)
1231   //
1232   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1233
1234   // integrate over regions
1235   // cache it for efficiency (usually more than one efficiency is requested)
1236
1237   AliCFContainer* sourceContainer = 0;
1238
1239   if (source == 0)
1240   {
1241     if (!fCache)
1242     {
1243       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1244       for (UInt_t i = 1; i < fkRegions; i++)
1245         if (fTrackHist[i])
1246           fCache->Add(fTrackHist[i]);
1247     }
1248     sourceContainer = fCache;
1249   }
1250   else if (source == 1 || source == 2)
1251   {
1252     sourceContainer = fTrackHistEfficiency;
1253     // step offset because we start with kCFStepAnaTopology
1254     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1255     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1256   }
1257   else
1258     return 0;
1259
1260   // reset all limits and set the right ones except those in axis1, axis2 and axis3
1261   ResetBinLimits(sourceContainer->GetGrid(step1));
1262   ResetBinLimits(sourceContainer->GetGrid(step2));
1263   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1264   {
1265     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1266     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1267     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1268   }
1269   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1270   {
1271     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1272     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1273   }
1274   if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1275   {
1276     sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1277     sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1278   }
1279
1280   TH1* measured = 0;
1281   TH1* generated = 0;
1282
1283   if (axis3 >= 0)
1284   {
1285     generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1286     measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1287   }
1288   else if (axis2 >= 0)
1289   {
1290     generated = sourceContainer->Project(step1, axis1, axis2);
1291     measured = sourceContainer->Project(step2, axis1, axis2);
1292   }
1293   else
1294   {
1295     generated = sourceContainer->Project(step1, axis1);
1296     measured = sourceContainer->Project(step2, axis1);
1297   }
1298
1299   // check for bins with less than 50 entries, print warning
1300   Int_t binBegin[3];
1301   Int_t binEnd[3];
1302
1303   binBegin[0] = 1;
1304   binBegin[1] = 1;
1305   binBegin[2] = 1;
1306
1307   binEnd[0] = generated->GetNbinsX();
1308   binEnd[1] = generated->GetNbinsY();
1309   binEnd[2] = generated->GetNbinsZ();
1310
1311   if (fEtaMax > fEtaMin)
1312   {
1313     if (axis1 == 0)
1314     {
1315       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1316       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
1317     }
1318     if (axis2 == 0)
1319     {
1320       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1321       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
1322     }
1323     if (axis3 == 0)
1324     {
1325       binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1326       binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
1327     }
1328   }
1329
1330   if (fPtMax > fPtMin)
1331   {
1332     // TODO this is just checking up to 15 for now
1333     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1334     if (axis1 == 1)
1335     {
1336       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1337       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
1338     }
1339     if (axis2 == 1)
1340     {
1341       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1342       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
1343     }
1344     if (axis3 == 1)
1345     {
1346       binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1347       binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
1348     }
1349   }
1350
1351   Int_t total = 0;
1352   Int_t count = 0;
1353   Int_t vars[3];
1354
1355   for (Int_t i=0; i<3; i++)
1356     vars[i] = binBegin[i];
1357
1358   const Int_t limit = 50;
1359   while (1)
1360   {
1361     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1362     {
1363       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1364       count++;
1365     }
1366     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1367     {
1368       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
1369         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1370         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1371         generated->GetBinContent(vars[0], vars[1]));
1372       count++;
1373     }
1374     else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1375     {
1376       Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
1377         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1378         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1379         generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1380         generated->GetBinContent(vars[0], vars[1], vars[2]));
1381       count++;
1382     }
1383
1384     vars[2]++;
1385     if (vars[2] == binEnd[2]+1)
1386     {
1387       vars[2] = binBegin[2];
1388       vars[1]++;
1389     }
1390
1391     if (vars[1] == binEnd[1]+1)
1392     {
1393       vars[1] = binBegin[1];
1394       vars[0]++;
1395     }
1396
1397     if (vars[0] == binEnd[0]+1)
1398       break;
1399     total++;
1400   }
1401
1402   Printf("Correction has %d empty bins (out of %d bins)", count, total);
1403
1404   // rebin if required
1405   if (source == 2)
1406   {
1407     TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1408
1409     if (axis->GetNbins() < measured->GetNbinsX())
1410     {
1411       if (axis2 != -1)
1412       {
1413         // 2d rebin with variable axis does not exist in root
1414
1415         TH1* tmp = measured;
1416         measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1417         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1418           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1419           {
1420             ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1421             measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1422           }
1423         delete tmp;
1424
1425         tmp = generated;
1426         generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1427         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1428           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1429           {
1430             ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1431             generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1432           }
1433         delete tmp;
1434       }
1435       else
1436       {
1437         TH1* tmp = measured;
1438         measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
1439         delete tmp;
1440
1441         tmp = generated;
1442         generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
1443         delete tmp;
1444       }
1445     }
1446     else if (axis->GetNbins() > measured->GetNbinsX())
1447     {
1448       if (axis2 != -1)
1449         AliFatal("Rebinning only works for 1d at present");
1450
1451       // this is an unfortunate case where the number of bins has to be increased in principle
1452       // there is a region where the binning is finner in one histogram and a region where it is the other way round
1453       // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
1454       // only a certain binning is supported here
1455       if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
1456         AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
1457
1458       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};
1459
1460       // reduce binning below 5 GeV/c
1461       TH1* tmp = measured;
1462       measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
1463       delete tmp;
1464
1465       // increase binning above 5 GeV/c
1466       tmp = measured;
1467       measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
1468       for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
1469       {
1470         measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
1471         measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
1472       }
1473       delete tmp;
1474
1475       // reduce binning below 5 GeV/c
1476       tmp = generated;
1477       generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
1478       delete tmp;
1479
1480       // increase binning above 5 GeV/c
1481       tmp = generated;
1482       generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
1483       for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
1484       {
1485         generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
1486         generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
1487       }
1488       delete tmp;
1489     }
1490   }
1491
1492   measured->Divide(measured, generated, 1, 1, "B");
1493
1494   delete generated;
1495
1496   ResetBinLimits(sourceContainer->GetGrid(step1));
1497   ResetBinLimits(sourceContainer->GetGrid(step2));
1498
1499   return measured;
1500 }
1501
1502 //____________________________________________________________________
1503 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
1504 {
1505   // creates a event-level efficiency by dividing step2 by step1
1506   // projected to axis1 and axis2 (optional if >= 0)
1507
1508   if (ptLeadMax > ptLeadMin)
1509   {
1510     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
1511     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
1512   }
1513
1514   TH1* measured = 0;
1515   TH1* generated = 0;
1516
1517   if (axis2 >= 0)
1518   {
1519     generated = fEventHist->Project(step1, axis1, axis2);
1520     measured = fEventHist->Project(step2, axis1, axis2);
1521   }
1522   else
1523   {
1524     generated = fEventHist->Project(step1, axis1);
1525     measured = fEventHist->Project(step2, axis1);
1526   }
1527
1528   measured->Divide(measured, generated, 1, 1, "B");
1529
1530   delete generated;
1531
1532   if (ptLeadMax > ptLeadMin)
1533   {
1534     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
1535     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
1536   }
1537
1538   return measured;
1539 }
1540
1541 //____________________________________________________________________
1542 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
1543 {
1544   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
1545   // where the matching is done of the z axis of hist1 with the x axis of hist2
1546
1547   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
1548     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
1549
1550   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
1551   {
1552     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
1553     {
1554       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
1555       {
1556         if (hist2->GetBinContent(z) > 0)
1557         {
1558           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
1559           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
1560         }
1561         else
1562         {
1563           hist1->SetBinContent(x, y, z, 0);
1564           hist1->SetBinError(x, y, z, 0);
1565         }
1566       }
1567     }
1568   }
1569 }
1570
1571 //____________________________________________________________________
1572 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
1573 {
1574   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
1575   // in the given region (sum over all regions is calculated if region == -1)
1576   // done by weighting the track-level distribution with the number of events as function of leading pT
1577   // and then calculating the ratio between the distributions
1578   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
1579   //   no projection is done if axis == 0
1580   // weighting: 0 = tracks weighted with events (as discussed above)
1581   //            1 = only track bias is returned
1582   //            2 = only event bias is returned
1583
1584   AliCFContainer* tmp = 0;
1585
1586   if (region == -1)
1587   {
1588     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
1589     for (UInt_t i = 1; i < fkRegions; i++)
1590       if (fTrackHist[i])
1591         tmp->Add(fTrackHist[i]);
1592   }
1593   else if (region == kMin && fCombineMinMax)
1594   {
1595     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
1596     tmp->Add(fTrackHist[kMax]);
1597   }
1598   else
1599     tmp = fTrackHist[region];
1600
1601   ResetBinLimits(tmp->GetGrid(step1));
1602   ResetBinLimits(fEventHist->GetGrid(step1));
1603   SetBinLimits(tmp->GetGrid(step1));
1604
1605   ResetBinLimits(tmp->GetGrid(step2));
1606   ResetBinLimits(fEventHist->GetGrid(step2));
1607   SetBinLimits(tmp->GetGrid(step2));
1608
1609   TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
1610   TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
1611   if (weighting == 0)
1612     WeightHistogram(hist1, events1);
1613
1614   TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
1615   TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
1616   if (weighting == 0)
1617     WeightHistogram(hist2, events2);
1618
1619   TH1* generated = hist1;
1620   TH1* measured = hist2;
1621
1622   if (weighting == 0 || weighting == 1)
1623   {
1624     if (axis)
1625     {
1626       if (leadPtMax > leadPtMin)
1627       {
1628         hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1629         hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
1630       }
1631
1632       if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
1633       {
1634         hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1635         hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
1636       }
1637
1638       generated = hist1->Project3D(axis);
1639       measured  = hist2->Project3D(axis);
1640
1641       // delete hists here if projection has been done
1642       delete hist1;
1643       delete hist2;
1644     }
1645     delete events1;
1646     delete events2;
1647   }
1648   else if (weighting == 2)
1649   {
1650     delete hist1;
1651     delete hist2;
1652     generated = events1;
1653     measured = events2;
1654   }
1655
1656   measured->Divide(generated);
1657
1658   delete generated;
1659
1660   ResetBinLimits(tmp->GetGrid(step1));
1661   ResetBinLimits(tmp->GetGrid(step2));
1662
1663   if ((region == -1) || (region == kMin && fCombineMinMax))
1664     delete tmp;
1665
1666   return measured;
1667 }
1668
1669 //____________________________________________________________________
1670 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
1671 {
1672   // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
1673
1674   if (!fTrackHist[region])
1675     return;
1676
1677   THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
1678
1679   Int_t var1 = 1;
1680   Int_t var2 = 2;
1681
1682   if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
1683     AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
1684
1685   if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
1686     AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
1687
1688   // optimized implementation
1689   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1690   {
1691     Int_t bins[6];
1692
1693     Double_t value = grid->GetBinContent(binIdx, bins);
1694     Double_t error = grid->GetBinError(binIdx);
1695
1696     // check eta and phi axes
1697     if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
1698       continue;
1699     if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
1700       continue;
1701
1702     value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
1703     error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
1704
1705     grid->SetBinContent(bins, value);
1706     grid->SetBinError(bins, error);
1707   }
1708
1709   Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
1710 }
1711
1712 //____________________________________________________________________
1713 TH2* AliUEHist::GetCorrelatedContamination()
1714 {
1715   // 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!)
1716
1717   Int_t step1 = kCFStepTrackedOnlyPrim;
1718   Int_t step2 = kCFStepTracked;
1719
1720   fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
1721   fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
1722   TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
1723
1724   fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
1725   fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
1726   TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
1727
1728   tracksStep1->Divide(tracksStep2);
1729
1730   TH1* triggersStep1 = fEventHist->Project(step1, 0);
1731   TH1* triggersStep2 = fEventHist->Project(step2, 0);
1732
1733   TH1* singleParticle = GetTrackingContamination(1);
1734
1735   for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
1736     for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
1737       if (singleParticle->GetBinContent(x) > 0 && triggersStep1->GetBinContent(y) > 0)
1738         tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
1739       else
1740         tracksStep1->SetBinContent(x, y, 0);
1741
1742   delete singleParticle;
1743   delete tracksStep2;
1744   delete triggersStep1;
1745   delete triggersStep2;
1746
1747   return tracksStep1;
1748 }
1749
1750 //____________________________________________________________________
1751 TH2D* AliUEHist::GetTrackingEfficiency()
1752 {
1753   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1754   // integrates over the regions and all other variables than pT and eta to increase the statistics
1755   //
1756   // returned histogram has to be deleted by the user
1757
1758   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
1759 }
1760
1761 //____________________________________________________________________
1762 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
1763 {
1764   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1765   // integrates over the regions and all other variables than pT, centrality to increase the statistics
1766   //
1767   // returned histogram has to be deleted by the user
1768
1769   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
1770 }
1771
1772 //____________________________________________________________________
1773 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
1774 {
1775   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
1776   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1777
1778   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
1779 }
1780
1781 //____________________________________________________________________
1782 TH2D* AliUEHist::GetTrackingCorrection()
1783 {
1784   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1785   // integrates over the regions and all other variables than pT and eta to increase the statistics
1786   //
1787   // returned histogram has to be deleted by the user
1788
1789   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
1790 }
1791
1792 //____________________________________________________________________
1793 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
1794 {
1795   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1796   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1797
1798   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
1799 }
1800
1801 //____________________________________________________________________
1802 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
1803 {
1804   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1805   // integrates over the regions and all other variables than pT and eta to increase the statistics
1806   //
1807   // returned histogram has to be deleted by the user
1808
1809   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
1810 }
1811
1812 //____________________________________________________________________
1813 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
1814 {
1815   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1816   // integrates over the regions and all other variables than pT and centrality to increase the statistics
1817   //
1818   // returned histogram has to be deleted by the user
1819
1820   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
1821 }
1822
1823 //____________________________________________________________________
1824 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
1825 {
1826   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1827   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1828
1829   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
1830 }
1831
1832 //____________________________________________________________________
1833 TH2D* AliUEHist::GetTrackingContamination()
1834 {
1835   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1836   // integrates over the regions and all other variables than pT and eta to increase the statistics
1837   //
1838   // returned histogram has to be deleted by the user
1839
1840   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
1841 }
1842
1843 //____________________________________________________________________
1844 TH2D* AliUEHist::GetTrackingContaminationCentrality()
1845 {
1846   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1847   // integrates over the regions and all other variables than pT and centrality to increase the statistics
1848   //
1849   // returned histogram has to be deleted by the user
1850
1851   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
1852 }
1853
1854 //____________________________________________________________________
1855 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
1856 {
1857   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1858   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1859
1860   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
1861 }
1862
1863 //____________________________________________________________________
1864 const char* AliUEHist::GetRegionTitle(Region region)
1865 {
1866   // returns the name of the given region
1867
1868   switch (region)
1869   {
1870     case kToward:
1871       return "Towards";
1872     case kAway:
1873       return "Away";
1874     case kMin:
1875       return (fCombineMinMax) ? "Transverse" : "Min";
1876     case kMax:
1877       return "Max";
1878   }
1879
1880   return 0;
1881 }
1882
1883 //____________________________________________________________________
1884 const char* AliUEHist::GetStepTitle(CFStep step)
1885 {
1886   // returns the name of the given step
1887
1888   switch (step)
1889   {
1890     case kCFStepAll:
1891       return "All events";
1892     case kCFStepTriggered:
1893       return "Triggered";
1894     case kCFStepVertex:
1895       return "Primary Vertex";
1896     case kCFStepAnaTopology:
1897       return "Required analysis topology";
1898     case kCFStepTrackedOnlyPrim:
1899       return "Tracked (matched MC, only primaries)";
1900     case kCFStepTracked:
1901       return "Tracked (matched MC, all)";
1902     case kCFStepReconstructed:
1903       return "Reconstructed";
1904     case kCFStepRealLeading:
1905       return "Correct leading particle identified";
1906     case kCFStepBiasStudy:
1907       return "Bias study applying tracking efficiency";
1908     case kCFStepBiasStudy2:
1909       return "Bias study applying tracking efficiency in two steps";
1910   }
1911
1912   return 0;
1913 }
1914
1915 //____________________________________________________________________
1916 void AliUEHist::CopyReconstructedData(AliUEHist* from)
1917 {
1918   // copies those histograms extracted from ESD to this object
1919
1920   // TODO at present only the pointers are copied
1921
1922   for (Int_t region=0; region<4; region++)
1923   {
1924     if (!fTrackHist[region])
1925       continue;
1926
1927     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
1928     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1929     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
1930   }
1931
1932   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
1933   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
1934   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
1935 }
1936
1937 //____________________________________________________________________
1938 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
1939 {
1940   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
1941
1942   Float_t fitRangeBegin = 5.01;
1943   Float_t fitRangeEnd = 14.99;
1944   Float_t extendRangeBegin = 10.01;
1945
1946   if (fTrackHistEfficiency->GetNVar() == 3)
1947   {
1948     TH1* obj = GetTrackingEfficiency(1);
1949
1950     if (verbose)
1951     {
1952       new TCanvas;
1953       obj->Draw();
1954     }
1955
1956     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
1957
1958     Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
1959
1960     obj = GetTrackingContamination(1);
1961
1962     if (verbose)
1963     {
1964       new TCanvas;
1965       obj->Draw();
1966     }
1967
1968     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
1969
1970     Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
1971
1972     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);
1973
1974     // extend for full pT range
1975     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
1976       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
1977         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
1978         {
1979
1980           Int_t bins[3];
1981           bins[0] = x;
1982           bins[1] = y;
1983           bins[2] = z;
1984
1985           fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
1986           fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
1987           fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
1988         }
1989   }
1990   else if (fTrackHistEfficiency->GetNVar() == 4)
1991   {
1992     // fit in centrality intervals of 20% for efficiency, one bin for contamination
1993     Float_t* trackingEff = 0;
1994     Float_t* trackingCont = 0;
1995     Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
1996     Int_t nCentralityBins = 5;
1997
1998     Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
1999
2000     // 0 = eff; 1 = cont
2001     for (Int_t caseNo = 0; caseNo < 2; caseNo++)
2002     {
2003       Float_t* target = 0;
2004       Int_t centralityBinsLocal = nCentralityBins;
2005
2006       if (caseNo == 0)
2007       {
2008         trackingEff = new Float_t[centralityBinsLocal];
2009         target = trackingEff;
2010       }
2011       else
2012       {
2013         centralityBinsLocal = 1;
2014         trackingCont = new Float_t[centralityBinsLocal];
2015         target = trackingCont;
2016       }
2017
2018       for (Int_t i=0; i<centralityBinsLocal; i++)
2019       {
2020         if (centralityBinsLocal == 1)
2021           SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
2022         else
2023           SetCentralityRange(centralityBins[i] + 0.1, centralityBins[i+1] - 0.1);
2024         TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
2025         if (verbose)
2026         {
2027           new TCanvas;
2028           proj->DrawCopy();
2029         }
2030         if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
2031           target[i] = proj->GetFunction("pol0")->GetParameter(0);
2032         else
2033           target[i] = 0;
2034         Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
2035       }
2036     }
2037
2038     // extend for full pT range
2039     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2040       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2041         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2042         {
2043           for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
2044           {
2045
2046             Int_t bins[4];
2047             bins[0] = x;
2048             bins[1] = y;
2049             bins[2] = z;
2050             bins[3] = z2;
2051
2052             Int_t z2Bin = 0;
2053             while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
2054               z2Bin++;
2055
2056             //Printf("%d %d", z2, z2Bin);
2057
2058             fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2059             fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
2060             if (trackingCont[0] > 0)
2061               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
2062             else
2063               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
2064           }
2065        }
2066
2067      delete[] trackingEff;
2068      delete[] trackingCont;
2069    }
2070
2071    SetCentralityRange(0, 0);
2072 }
2073
2074 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
2075 {
2076   // corrects the dphi distribution with an extra factor close to dphi ~ 0
2077
2078   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
2079
2080   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
2081
2082   // optimized implementation
2083   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2084   {
2085     Int_t bins[6];
2086     Double_t value = grid->GetBinContent(binIdx, bins);
2087     Double_t error = grid->GetBinError(binIdx);
2088
2089     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
2090     if (TMath::Abs(binCenter) < 0.2)
2091     {
2092       value *= 0.985;
2093       error *= 0.985;
2094     }
2095     else if (TMath::Abs(binCenter) < 0.3)
2096     {
2097       value *= 0.9925;
2098       error *= 0.9925;
2099     }
2100
2101     grid->SetBinContent(bins, value);
2102     grid->SetBinError(bins, error);
2103   }
2104 }
2105
2106 void AliUEHist::Scale(Double_t factor)
2107 {
2108   // scales all contained histograms by the given factor
2109
2110   for (Int_t i=0; i<4; i++)
2111     if (fTrackHist[i])
2112       fTrackHist[i]->Scale(factor);
2113
2114   fEventHist->Scale(factor);
2115   fTrackHistEfficiency->Scale(factor);
2116 }
2117
2118 void AliUEHist::Reset()
2119 {
2120   // resets all contained histograms
2121
2122   for (Int_t i=0; i<4; i++)
2123     if (fTrackHist[i])
2124       for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
2125         fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2126
2127   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2128     fEventHist->GetGrid(step)->GetGrid()->Reset();
2129
2130   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
2131     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
2132 }