e08f9f1367eabb10159b61271c626f29dfa20e6c
[u/mrichter/AliRoot.git] / 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 }