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
  38 ClassImp(AliUEHist)
  39
  40 const Int_t AliUEHist::fgkCFSteps = 10;
  41
  42 AliUEHist::AliUEHist(const char* reqHist) :
  43   TObject(),
  44   fkRegions(4),
  45   fEventHist(0),
  46   fEtaMin(0),
  47   fEtaMax(0),
  48   fPtMin(0),
  49   fPtMax(0),
  50   fContaminationEnhancement(0),
  51   fCombineMinMax(0),
  52   fCache(0)
  53 {
  54   // Constructor
  55
  56   for (Int_t i=0; i<fkRegions; i++)
  57     fTrackHist[i] = 0;
  58
  59   if (strlen(reqHist) == 0)
  60     return;
  61
  62   const char* title = "";
  63
  64   // track level
  65   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
  66   Int_t iTrackBin[5];
  67   Double_t* trackBins[5];
  68   const char* trackAxisTitle[5];
  69
  70   // eta
  71   iTrackBin[0] = 20;
  72   Double_t etaBins[20+1];
  73   for (Int_t i=0; i<=iTrackBin[0]; i++)
  74     etaBins[i] = -1.0 + 0.1 * i;
  75   trackBins[0] = etaBins;
  76   trackAxisTitle[0] = "#eta";
  77
  78   // pT
  79   iTrackBin[1] = 39;
  80   Double_t pTBins[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
  81   trackBins[1] = pTBins;
  82   trackAxisTitle[1] = "p_{T} (GeV/c)";
  83
  84   // pT,lead binning 1
  85   const Int_t kNLeadingpTBins = 100;
  86   Double_t leadingpTBins[kNLeadingpTBins+1];
  87   for (Int_t i=0; i<=kNLeadingpTBins; i++)
  88     leadingpTBins[i] = 0.5 * i;
  89
  90   // pT,lead binning 2
  91   const Int_t kNLeadingpTBins2 = 10;
  92   Double_t leadingpTBins2[kNLeadingpTBins2+1];
  93   for (Int_t i=0; i<=kNLeadingpTBins2; i++)
  94     leadingpTBins2[i] = 5.0 * i;
  95
  96   // phi,lead
  97   const Int_t kNLeadingPhiBins = 40;
  98   Double_t leadingPhiBins[kNLeadingPhiBins+1];
  99   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
 100     leadingPhiBins[i] = -1.5 * TMath::Pi() + 1.0 / 40 * i * TMath::TwoPi();
 101
 102   // multiplicity
 103   const Int_t kNMultiplicityBins = 15;
 104   Double_t multiplicityBins[kNMultiplicityBins+1];
 105   for (Int_t i=0; i<=kNMultiplicityBins; i++)
 106     multiplicityBins[i] = -0.5 + i;
 107   multiplicityBins[kNMultiplicityBins] = 200;
 108
 109   trackBins[3] = multiplicityBins;
 110   iTrackBin[3] = kNMultiplicityBins;
 111   trackAxisTitle[3] = "multiplicity";
 112
 113   // selection depending on requested histogram
 114   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
 115   if (strcmp(reqHist, "NumberDensitypT") == 0)
 116   {
 117     axis = 0;
 118     title = "d^{2}N_{ch}/d#phid#eta";
 119   }
 120   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
 121   {
 122     axis = 1;
 123     title = "d^{2}N_{ch}/d#phid#eta";
 124   }
 125   else if (strcmp(reqHist, "SumpT") == 0)
 126   {
 127     axis = 0;
 128     title = "d^{2}#Sigma p_{T}/d#phid#eta";
 129   }
 130   else
 131     AliFatal(Form("Invalid histogram requested: %s", reqHist));
 132
 133   Int_t initRegions = fkRegions;
 134
 135   if (axis == 0)
 136   {
 137     trackBins[2] = leadingpTBins;
 138     iTrackBin[2] = kNLeadingpTBins;
 139     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 140
 141   }
 142   else if (axis == 1)
 143   {
 144     nTrackVars = 5;
 145     initRegions = 1;
 146
 147     iTrackBin[2] = kNLeadingpTBins2;
 148     trackBins[2] = leadingpTBins2;
 149     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
 150
 151     iTrackBin[4] = kNLeadingPhiBins;
 152     trackBins[4] = leadingPhiBins;
 153     trackAxisTitle[4] = "#phi w.r.t leading track";
 154   }
 155
 156   for (Int_t i=0; i<initRegions; i++)
 157   {
 158     fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
 159
 160     for (Int_t j=0; j<nTrackVars; j++)
 161     {
 162       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
 163       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
 164     }
 165
 166     SetStepNames(fTrackHist[i]);
 167   }
 168
 169   // track 3rd and 4th axis --> event 1st and 2nd axis
 170   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
 171
 172   fEventHist->SetBinLimits(0, trackBins[2]);
 173   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
 174
 175   fEventHist->SetBinLimits(1, trackBins[3]);
 176   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
 177
 178   SetStepNames(fEventHist);
 179 }
 180
 181 //_____________________________________________________________________________
 182 AliUEHist::AliUEHist(const AliUEHist &c) :
 183   TObject(),
 184   fkRegions(4),
 185   fEventHist(0),
 186   fEtaMin(0),
 187   fEtaMax(0),
 188   fPtMin(0),
 189   fPtMax(0),
 190   fContaminationEnhancement(0),
 191   fCombineMinMax(0),
 192   fCache(0)
 193 {
 194   //
 195   // AliUEHist copy constructor
 196   //
 197
 198   ((AliUEHist &) c).Copy(*this);
 199 }
 200
 201 //____________________________________________________________________
 202 void AliUEHist::SetStepNames(AliCFContainer* container)
 203 {
 204   // sets the names of the correction steps
 205
 206   for (Int_t i=0; i<fgkCFSteps; i++)
 207     container->SetStepTitle(i, GetStepTitle((CFStep) i));
 208 }
 209
 210 //____________________________________________________________________
 211 AliUEHist::~AliUEHist()
 212 {
 213   // Destructor
 214
 215   for (Int_t i=0; i<fkRegions; i++)
 216   {
 217     if (fTrackHist[i])
 218     {
 219       delete fTrackHist[i];
 220       fTrackHist[i] = 0;
 221     }
 222   }
 223
 224   if (fEventHist)
 225   {
 226     delete fEventHist;
 227     fEventHist = 0;
 228   }
 229
 230   if (fCache)
 231   {
 232     delete fCache;
 233     fCache = 0;
 234   }
 235 }
 236
 237 //____________________________________________________________________
 238 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
 239 {
 240   // assigment operator
 241
 242   if (this != &c)
 243     ((AliUEHist &) c).Copy(*this);
 244
 245   return *this;
 246 }
 247
 248 //____________________________________________________________________
 249 void AliUEHist::Copy(TObject& c) const
 250 {
 251   // copy function
 252
 253   AliUEHist& target = (AliUEHist &) c;
 254
 255   for (Int_t i=0; i<fkRegions; i++)
 256     if (fTrackHist[i])
 257       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
 258
 259   if (fEventHist)
 260     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
 261 }
 262
 263 //____________________________________________________________________
 264 Long64_t AliUEHist::Merge(TCollection* list)
 265 {
 266   // Merge a list of AliUEHist objects with this (needed for
 267   // PROOF).
 268   // Returns the number of merged objects (including this).
 269
 270   if (!list)
 271     return 0;
 272
 273   if (list->IsEmpty())
 274     return 1;
 275
 276   TIterator* iter = list->MakeIterator();
 277   TObject* obj;
 278
 279   // collections of objects
 280   const Int_t kMaxLists = fkRegions+1;
 281   TList** lists = new TList*[kMaxLists];
 282
 283   for (Int_t i=0; i<kMaxLists; i++)
 284     lists[i] = new TList;
 285
 286   Int_t count = 0;
 287   while ((obj = iter->Next())) {
 288
 289     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
 290     if (entry == 0)
 291       continue;
 292
 293     for (Int_t i=0; i<fkRegions; i++)
 294       if (entry->fTrackHist[i])
 295         lists[i]->Add(entry->fTrackHist[i]);
 296
 297     lists[fkRegions]->Add(entry->fEventHist);
 298
 299     count++;
 300   }
 301   for (Int_t i=0; i<fkRegions; i++)
 302     if (fTrackHist[i])
 303       fTrackHist[i]->Merge(lists[i]);
 304
 305   fEventHist->Merge(lists[fkRegions]);
 306
 307   for (Int_t i=0; i<kMaxLists; i++)
 308     delete lists[i];
 309
 310   delete[] lists;
 311
 312   return count+1;
 313 }
 314
 315 //____________________________________________________________________
 316 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
 317 {
 318   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
 319
 320   if (fEtaMax > fEtaMin)
 321     grid->SetRangeUser(0, fEtaMin, fEtaMax);
 322   if (fPtMax > fPtMin)
 323     grid->SetRangeUser(1, fPtMin, fPtMax);
 324 }
 325
 326 //____________________________________________________________________
 327 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
 328 {
 329   // resets all bin limits
 330
 331   for (Int_t i=0; i<grid->GetNVar(); i++)
 332     if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange))
 333       grid->SetRangeUser(i, 0, -1);
 334 }
 335
 336 //____________________________________________________________________
 337 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
 338 {
 339   // prints the number of empty bins in the track end event histograms in the given step
 340
 341   Int_t binBegin[4];
 342   Int_t binEnd[4];
 343
 344   for (Int_t i=0; i<4; i++)
 345   {
 346     binBegin[i] = 1;
 347     binEnd[i]   = fTrackHist[0]->GetNBins(i);
 348   }
 349
 350   if (fEtaMax > fEtaMin)
 351   {
 352     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
 353     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
 354   }
 355
 356   if (fPtMax > fPtMin)
 357   {
 358     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
 359     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
 360   }
 361
 362   if (ptLeadMax > ptLeadMin)
 363   {
 364     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
 365     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
 366   }
 367
 368   // start from multiplicity 1
 369   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
 370
 371   for (Int_t region=0; region<fkRegions; region++)
 372   {
 373     Int_t total = 0;
 374     Int_t count = 0;
 375     Int_t vars[4];
 376
 377     for (Int_t i=0; i<4; i++)
 378       vars[i] = binBegin[i];
 379
 380     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
 381     while (1)
 382     {
 383       if (grid->GetElement(vars) == 0)
 384       {
 385         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
 386           grid->GetBinCenter(0, vars[0]),
 387           grid->GetBinCenter(1, vars[1]),
 388           grid->GetBinCenter(2, vars[2]),
 389           grid->GetBinCenter(3, vars[3])
 390         );
 391         count++;
 392       }
 393
 394       vars[3]++;
 395       for (Int_t i=3; i>0; i--)
 396       {
 397         if (vars[i] == binEnd[i]+1)
 398         {
 399           vars[i] = binBegin[i];
 400           vars[i-1]++;
 401         }
 402       }
 403
 404       if (vars[0] == binEnd[0]+1)
 405         break;
 406       total++;
 407     }
 408
 409     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
 410   }
 411 }
 412
 413 //____________________________________________________________________
 414 TH1D* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax)
 415 {
 416   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
 417   //
 418   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
 419   // Histogram has to be deleted by the caller of the function
 420
 421   // unzoom all axes
 422   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 423   ResetBinLimits(fEventHist->GetGrid(step));
 424
 425   SetBinLimits(fTrackHist[region]->GetGrid(step));
 426
 427   TH1D* tracks = 0;
 428
 429   if (ptLeadMin < 0)
 430   {
 431     tracks = fTrackHist[region]->ShowProjection(2, step);
 432     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
 433     if (fCombineMinMax && region == kMin)
 434     {
 435       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 436       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
 437
 438       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
 439       tracks->Add(tracks2);
 440
 441       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
 442     }
 443
 444     // normalize to get a density (deta dphi)
 445     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 446     Float_t phiRegion = TMath::TwoPi() / 3;
 447     if (!fCombineMinMax && region == kMin)
 448       phiRegion /= 2;
 449     Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 450     //Printf("Normalization: %f", normalization);
 451     tracks->Scale(1.0 / normalization);
 452
 453     TH1D* events = fEventHist->ShowProjection(0, step);
 454     tracks->Divide(events);
 455   }
 456   else
 457   {
 458     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, ptLeadMin, ptLeadMax);
 459     tracks = fTrackHist[region]->GetGrid(step)->Project(4);
 460     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
 461
 462     // normalize to get a density (deta dphi)
 463     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
 464     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1) * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
 465     //Printf("Normalization: %f", normalization);
 466     tracks->Scale(1.0 / normalization);
 467
 468     TH1D* events = fEventHist->ShowProjection(0, step);
 469     Int_t nEvents = (Int_t) events->Integral(events->FindBin(ptLeadMin), events->FindBin(ptLeadMax));
 470     if (nEvents > 0)
 471       tracks->Scale(1.0 / nEvents);
 472   }
 473
 474   ResetBinLimits(fTrackHist[region]->GetGrid(step));
 475
 476   return tracks;
 477 }
 478
 479 //____________________________________________________________________
 480 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
 481 {
 482   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
 483   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
 484   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
 485   //
 486   // if trackCorrection is 0, just copies content from step1 to step2
 487
 488   for (Int_t region=0; region<fkRegions; region++)
 489     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
 490 }
 491
 492 //____________________________________________________________________
 493 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
 494 {
 495   //
 496   // see documentation of CorrectTracks above
 497   //
 498
 499   if (!fTrackHist[region])
 500     return;
 501
 502   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
 503   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
 504
 505   // clear target histogram
 506   target->Reset();
 507
 508   if (trackCorrection != 0)
 509   {
 510     if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
 511       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
 512
 513     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
 514       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
 515   }
 516
 517   // optimized implementation
 518   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
 519   {
 520     Int_t bins[5];
 521     Double_t value = grid->GetBinContent(binIdx, bins);
 522     Double_t error = grid->GetBinError(binIdx);
 523
 524     if (trackCorrection != 0)
 525     {
 526       if (var2 < 0)
 527       {
 528         value *= trackCorrection->GetBinContent(bins[var1]);
 529         error *= trackCorrection->GetBinContent(bins[var1]);
 530       }
 531       else
 532       {
 533         value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 534         error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
 535       }
 536     }
 537
 538     target->SetBinContent(bins, value);
 539     target->SetBinError(bins, error);
 540   }
 541 }
 542
 543 //____________________________________________________________________
 544 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection, Int_t var)
 545 {
 546   // corrects from step1 to step2 by multiplying the events with eventCorrection
 547   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
 548   //
 549   // if eventCorrection is 0, just copies content from step1 to step2
 550
 551   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
 552   AliCFGridSparse* target = fEventHist->GetGrid(step2);
 553
 554   // clear target histogram
 555   target->GetGrid()->Reset();
 556
 557   if (eventCorrection != 0 && grid->GetNBins(var) != eventCorrection->GetNbinsX())
 558     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var), eventCorrection->GetNbinsX()));
 559
 560   Int_t bins[2];
 561   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
 562   {
 563     bins[0] = x;
 564     for (Int_t y = 1; y <= grid->GetNBins(1); y++)
 565     {
 566       bins[1] = y;
 567
 568       Double_t value = grid->GetElement(bins);
 569       if (value != 0)
 570       {
 571         Double_t error = grid->GetElementError(bins);
 572
 573         if (eventCorrection != 0)
 574         {
 575           value *= eventCorrection->GetBinContent(bins[var]);
 576           error *= eventCorrection->GetBinContent(bins[var]);
 577         }
 578
 579         target->SetElement(bins, value);
 580         target->SetElementError(bins, error);
 581       }
 582     }
 583   }
 584 }
 585
 586 //____________________________________________________________________
 587 void AliUEHist::Correct(AliUEHist* corrections)
 588 {
 589   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
 590   //
 591   // in this object the data is expected in the step kCFStepReconstructed
 592
 593   // ---- track level
 594
 595   // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
 596   // extracted as function of leading pT
 597   for (Int_t region = 0; region < fkRegions; region++)
 598   {
 599     if (!fTrackHist[region])
 600       continue;
 601
 602     //TH1D* leadingBias = (TH1D*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, "z"); // from MC
 603     TH1D* leadingBias = (TH1D*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, "z");          // from data
 604     CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBias, 2);
 605     if (region == kMin && fCombineMinMax)
 606     {
 607       CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBias, 2);
 608       delete leadingBias;
 609       break;
 610     }
 611     delete leadingBias;
 612   }
 613   CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, 0);
 614
 615   // correct with kCFStepTracked --> kCFStepTrackedOnlyPrim
 616   TH2D* contamination = corrections->GetTrackingContamination();
 617   if (corrections->fContaminationEnhancement)
 618   {
 619     Printf("Applying contamination enhancement");
 620
 621     for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
 622       for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
 623         contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
 624   }
 625   CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
 626   CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
 627   delete contamination;
 628
 629   // correct with kCFStepTrackedOnlyPrim --> kCFStepAnaTopology
 630   TH2D* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
 631   CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0, 1);
 632   CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 0);
 633   delete efficiencyCorrection;
 634
 635   // copy
 636   CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
 637   CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
 638
 639   // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
 640   // practically independent of low pT cut
 641   TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
 642
 643   // convert stage from true multiplicity to observed multiplicity by simple conversion factor
 644   TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
 645   vertexCorrectionObs->Reset();
 646
 647   TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
 648   vertexCorrection->Fit(func, "0", "", 0.8, 4);
 649
 650   for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
 651   {
 652     Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
 653     if (xPos < 4)
 654       vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
 655     else
 656       vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
 657   }
 658
 659   new TCanvas;
 660   vertexCorrection->DrawCopy();
 661   vertexCorrectionObs->SetLineColor(2);
 662   vertexCorrectionObs->DrawCopy("same");
 663
 664   CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
 665   CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
 666   delete vertexCorrectionObs;
 667   delete vertexCorrection;
 668
 669   // copy
 670   CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
 671   CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
 672 }
 673
 674 //____________________________________________________________________
 675 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2)
 676 {
 677   // creates a track-level efficiency by dividing step2 by step1
 678   // projected to axis1 and axis2 (optional if >= 0)
 679
 680   // integrate over regions
 681   // cache it for efficiency (usually more than one efficiency is requested)
 682   if (!fCache)
 683   {
 684     fCache = (AliCFContainer*) fTrackHist[0]->Clone();
 685     for (Int_t i = 1; i < fkRegions; i++)
 686       if (fTrackHist[i])
 687         fCache->Add(fTrackHist[i]);
 688   }
 689
 690   // reset all limits and set the right ones except those in axis1 and axis2
 691   ResetBinLimits(fCache->GetGrid(step1));
 692   ResetBinLimits(fCache->GetGrid(step2));
 693   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0)
 694   {
 695     fCache->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
 696     fCache->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
 697   }
 698   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1)
 699   {
 700     fCache->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
 701     fCache->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
 702   }
 703
 704   TH1* measured = 0;
 705   TH1* generated = 0;
 706
 707   if (axis2 >= 0)
 708   {
 709     generated = fCache->Project(axis1, axis2, step1);
 710     measured = fCache->Project(axis1, axis2, step2);
 711   }
 712   else
 713   {
 714     generated = fCache->Project(axis1, step1);
 715     measured = fCache->Project(axis1, step2);
 716   }
 717
 718   // check for bins with less than 100 entries, print warning
 719   Int_t binBegin[2];
 720   Int_t binEnd[2];
 721
 722   binBegin[0] = 1;
 723   binBegin[1] = 1;
 724
 725   binEnd[0] = generated->GetNbinsX();
 726   binEnd[1] = generated->GetNbinsY();
 727
 728   if (fEtaMax > fEtaMin)
 729   {
 730     if (axis1 == 0)
 731     {
 732       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
 733       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
 734     }
 735     if (axis2 == 0)
 736     {
 737       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
 738       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
 739     }
 740   }
 741
 742   if (fPtMax > fPtMin)
 743   {
 744     // TODO this is just checking up to 15 for now
 745     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
 746     if (axis1 == 1)
 747     {
 748       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
 749       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
 750     }
 751     if (axis2 == 1)
 752     {
 753       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
 754       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
 755     }
 756   }
 757
 758   Int_t total = 0;
 759   Int_t count = 0;
 760   Int_t vars[2];
 761
 762   for (Int_t i=0; i<2; i++)
 763     vars[i] = binBegin[i];
 764
 765   const Int_t limit = 50;
 766   while (1)
 767   {
 768     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
 769     {
 770       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
 771       count++;
 772     }
 773     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
 774     {
 775       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
 776         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
 777         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
 778         generated->GetBinContent(vars[0], vars[1]));
 779       count++;
 780     }
 781
 782     vars[1]++;
 783     if (vars[1] == binEnd[1]+1)
 784     {
 785       vars[1] = binBegin[1];
 786       vars[0]++;
 787     }
 788
 789     if (vars[0] == binEnd[0]+1)
 790       break;
 791     total++;
 792   }
 793
 794   Printf("Correction has %d empty bins (out of %d bins)", count, total);
 795
 796   measured->Divide(measured, generated, 1, 1, "B");
 797
 798   delete generated;
 799
 800   return measured;
 801 }
 802
 803 //____________________________________________________________________
 804 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
 805 {
 806   // creates a event-level efficiency by dividing step2 by step1
 807   // projected to axis1 and axis2 (optional if >= 0)
 808
 809   if (ptLeadMax > ptLeadMin)
 810   {
 811     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
 812     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
 813   }
 814
 815   TH1* measured = 0;
 816   TH1* generated = 0;
 817
 818   if (axis2 >= 0)
 819   {
 820     generated = fEventHist->Project(axis1, axis2, step1);
 821     measured = fEventHist->Project(axis1, axis2, step2);
 822   }
 823   else
 824   {
 825     generated = fEventHist->Project(axis1, step1);
 826     measured = fEventHist->Project(axis1, step2);
 827   }
 828
 829   measured->Divide(measured, generated, 1, 1, "B");
 830
 831   delete generated;
 832
 833   if (ptLeadMax > ptLeadMin)
 834   {
 835     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
 836     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
 837   }
 838
 839   return measured;
 840 }
 841
 842 //____________________________________________________________________
 843 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
 844 {
 845   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
 846   // where the matching is done of the z axis of hist1 with the x axis of hist2
 847
 848   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
 849     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
 850
 851   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
 852   {
 853     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
 854     {
 855       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
 856       {
 857         if (hist2->GetBinContent(z) > 0)
 858         {
 859           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
 860           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
 861         }
 862         else
 863         {
 864           hist1->SetBinContent(x, y, z, 0);
 865           hist1->SetBinError(x, y, z, 0);
 866         }
 867       }
 868     }
 869   }
 870 }
 871
 872 //____________________________________________________________________
 873 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax)
 874 {
 875   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
 876   // in the given region (sum over all regions is calculated if region == -1)
 877   // done by weighting the track-level distribution with the number of events as function of leading pT
 878   // and then calculating the ratio between the distributions
 879   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
 880   //   no projection is done if axis == 0
 881
 882   AliCFContainer* tmp = 0;
 883
 884   if (region == -1)
 885   {
 886     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
 887     for (Int_t i = 1; i < fkRegions; i++)
 888       if (fTrackHist[i])
 889         tmp->Add(fTrackHist[i]);
 890   }
 891   else if (region == kMin && fCombineMinMax)
 892   {
 893     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
 894     tmp->Add(fTrackHist[kMax]);
 895   }
 896   else
 897     tmp = fTrackHist[region];
 898
 899   ResetBinLimits(tmp->GetGrid(step1));
 900   ResetBinLimits(fEventHist->GetGrid(step1));
 901   SetBinLimits(tmp->GetGrid(step1));
 902
 903   ResetBinLimits(tmp->GetGrid(step2));
 904   ResetBinLimits(fEventHist->GetGrid(step2));
 905   SetBinLimits(tmp->GetGrid(step2));
 906
 907   TH1D* events1 = fEventHist->Project(0, step1);
 908   TH3D* hist1 = tmp->Project(0, 1, 2, step1);
 909   WeightHistogram(hist1, events1);
 910
 911   TH1D* events2 = fEventHist->Project(0, step2);
 912   TH3D* hist2 = tmp->Project(0, 1, 2, step2);
 913   WeightHistogram(hist2, events2);
 914
 915   TH1* generated = hist1;
 916   TH1* measured = hist2;
 917
 918   if (axis)
 919   {
 920     if (leadPtMax > leadPtMin)
 921     {
 922       hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
 923       hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
 924     }
 925
 926     if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
 927     {
 928       hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
 929       hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
 930     }
 931     if (fPtMax > fPtMin && !TString(axis).Contains("y"))
 932     {
 933       hist1->GetYaxis()->SetRangeUser(fPtMin, fPtMax);
 934       hist2->GetYaxis()->SetRangeUser(fPtMin, fPtMax);
 935     }
 936
 937     generated = hist1->Project3D(axis);
 938     measured  = hist2->Project3D(axis);
 939
 940     // delete hists here if projection has been done
 941     delete hist1;
 942     delete hist2;
 943   }
 944
 945   measured->Divide(generated);
 946
 947   delete events1;
 948   delete events2;
 949   delete generated;
 950
 951   ResetBinLimits(tmp->GetGrid(step1));
 952   ResetBinLimits(tmp->GetGrid(step2));
 953
 954   if ((region == -1) || (region == kMin && fCombineMinMax))
 955     delete tmp;
 956
 957   return measured;
 958 }
 959
 960 //____________________________________________________________________
 961 TH2D* AliUEHist::GetTrackingEfficiency()
 962 {
 963   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
 964   // integrates over the regions and all other variables than pT and eta to increase the statistics
 965   //
 966   // returned histogram has to be deleted by the user
 967
 968   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
 969 }
 970
 971 //____________________________________________________________________
 972 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
 973 {
 974   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
 975   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
 976
 977   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
 978 }
 979
 980 //____________________________________________________________________
 981 TH2D* AliUEHist::GetTrackingCorrection()
 982 {
 983   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
 984   // integrates over the regions and all other variables than pT and eta to increase the statistics
 985   //
 986   // returned histogram has to be deleted by the user
 987
 988   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
 989 }
 990
 991 //____________________________________________________________________
 992 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
 993 {
 994   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
 995   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
 996
 997   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
 998 }
 999
1000 //____________________________________________________________________
1001 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
1002 {
1003   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1004   // integrates over the regions and all other variables than pT and eta to increase the statistics
1005   //
1006   // returned histogram has to be deleted by the user
1007
1008   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
1009 }
1010
1011 //____________________________________________________________________
1012 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
1013 {
1014   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
1015   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1016
1017   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
1018 }
1019
1020 //____________________________________________________________________
1021 TH2D* AliUEHist::GetTrackingContamination()
1022 {
1023   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1024   // integrates over the regions and all other variables than pT and eta to increase the statistics
1025   //
1026   // returned histogram has to be deleted by the user
1027
1028   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
1029 }
1030
1031 //____________________________________________________________________
1032 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
1033 {
1034   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
1035   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
1036
1037   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
1038 }
1039
1040 //____________________________________________________________________
1041 const char* AliUEHist::GetRegionTitle(Region region)
1042 {
1043   // returns the name of the given region
1044
1045   switch (region)
1046   {
1047     case kToward:
1048       return "Towards";
1049     case kAway:
1050       return "Away";
1051     case kMin:
1052       return (fCombineMinMax) ? "Transverse" : "Min";
1053     case kMax:
1054       return "Max";
1055   }
1056
1057   return 0;
1058 }
1059
1060 //____________________________________________________________________
1061 const char* AliUEHist::GetStepTitle(CFStep step)
1062 {
1063   // returns the name of the given step
1064
1065   switch (step)
1066   {
1067     case kCFStepAll:
1068       return "All events";
1069     case kCFStepTriggered:
1070       return "Triggered";
1071     case kCFStepVertex:
1072       return "Primary Vertex";
1073     case kCFStepAnaTopology:
1074       return "Required analysis topology";
1075     case kCFStepTrackedOnlyPrim:
1076       return "Tracked (matched MC, only primaries)";
1077     case kCFStepTracked:
1078       return "Tracked (matched MC, all)";
1079     case kCFStepReconstructed:
1080       return "Reconstructed";
1081     case kCFStepRealLeading:
1082       return "Correct leading particle identified";
1083     case kCFStepBiasStudy:
1084       return "Bias study applying tracking efficiency";
1085     case kCFStepBiasStudy2:
1086       return "Bias study applying tracking efficiency in two steps";
1087   }
1088
1089   return 0;
1090 }