From Sara and Jan Fiete:
[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
37 ClassImp(AliUEHist)
38
39 const Int_t AliUEHist::fgkCFSteps = 10;
40
41 AliUEHist::AliUEHist(const char* reqHist) : 
42   fkRegions(4),
43   fEventHist(0),
44   fEtaMin(0),
45   fEtaMax(0),
46   fPtMin(0),
47   fPtMax(0),
48   fCombineMinMax(0),
49   fCache(0)
50 {
51   // Constructor
52     
53   for (Int_t i=0; i<fkRegions; i++)
54     fTrackHist[i] = 0;
55     
56   if (strlen(reqHist) == 0)
57     return;
58     
59   const char* title = "";
60     
61   // track level
62   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
63   Int_t iTrackBin[5];
64   Double_t* trackBins[5];
65   const char* trackAxisTitle[5];
66   
67   // eta
68   iTrackBin[0] = 20;
69   Double_t etaBins[20+1];
70   for (Int_t i=0; i<=iTrackBin[0]; i++)
71     etaBins[i] = -1.0 + 0.1 * i;
72   trackBins[0] = etaBins;
73   trackAxisTitle[0] = "#eta";
74   
75   // pT
76   iTrackBin[1] = 39;
77   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};
78   trackBins[1] = pTBins;
79   trackAxisTitle[1] = "p_{T} (GeV/c)";
80   
81   // pT,lead binning 1
82   const Int_t kNLeadingpTBins = 100;
83   Double_t leadingpTBins[kNLeadingpTBins+1];
84   for (Int_t i=0; i<=kNLeadingpTBins; i++)
85     leadingpTBins[i] = 0.5 * i;
86   
87   // pT,lead binning 2
88   const Int_t kNLeadingpTBins2 = 10;
89   Double_t leadingpTBins2[kNLeadingpTBins2+1];
90   for (Int_t i=0; i<=kNLeadingpTBins2; i++)
91     leadingpTBins2[i] = 5.0 * i;
92   
93   // phi,lead
94   const Int_t kNLeadingPhiBins = 40;
95   Double_t leadingPhiBins[kNLeadingPhiBins+1];
96   for (Int_t i=0; i<=kNLeadingPhiBins; i++)
97     leadingPhiBins[i] = -1.5 * TMath::Pi() + 1.0 / 40 * i * TMath::TwoPi();
98     
99   // multiplicity
100   const Int_t kNMultiplicityBins = 15;
101   Double_t multiplicityBins[kNMultiplicityBins+1];
102   for (Int_t i=0; i<=kNMultiplicityBins; i++)
103     multiplicityBins[i] = -0.5 + i;
104   multiplicityBins[kNMultiplicityBins] = 200;
105   
106   trackBins[3] = multiplicityBins;
107   iTrackBin[3] = kNMultiplicityBins;
108   trackAxisTitle[3] = "multiplicity";
109   
110   // selection depending on requested histogram
111   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
112   if (strcmp(reqHist, "NumberDensitypT") == 0)
113   {
114     axis = 0;
115     title = "d^{2}N_{ch}/d#phid#eta";
116   }
117   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
118   {
119     axis = 1;
120     title = "d^{2}N_{ch}/d#phid#eta";
121   }
122   else if (strcmp(reqHist, "SumpT") == 0)
123   {
124     axis = 0;
125     title = "d^{2}#Sigma p_{T}/d#phid#eta";
126   }
127   else
128     AliFatal(Form("Invalid histogram requested: %s", reqHist));
129   
130   Int_t initRegions = fkRegions;
131   
132   if (axis == 0)
133   {
134     trackBins[2] = leadingpTBins;
135     iTrackBin[2] = kNLeadingpTBins;
136     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
137     
138   }
139   else if (axis == 1)
140   {
141     nTrackVars = 5;
142     initRegions = 1;
143   
144     iTrackBin[2] = kNLeadingpTBins2;
145     trackBins[2] = leadingpTBins2;
146     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
147     
148     iTrackBin[4] = kNLeadingPhiBins;
149     trackBins[4] = leadingPhiBins;
150     trackAxisTitle[4] = "#phi w.r.t leading track";
151   }
152     
153   for (Int_t i=0; i<initRegions; i++)
154   {
155     fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
156     
157     for (Int_t j=0; j<nTrackVars; j++)
158     {
159       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
160       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
161     }
162     
163     SetStepNames(fTrackHist[i]);
164   }
165   
166   // track 3rd and 4th axis --> event 1st and 2nd axis
167   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2);
168   
169   fEventHist->SetBinLimits(0, trackBins[2]);
170   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
171   
172   fEventHist->SetBinLimits(1, trackBins[3]);
173   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
174   
175   SetStepNames(fEventHist);
176 }
177
178 //____________________________________________________________________
179 void AliUEHist::SetStepNames(AliCFContainer* container)
180 {
181   // sets the names of the correction steps
182   
183   for (Int_t i=0; i<fgkCFSteps; i++)
184     container->SetStepTitle(i, GetStepTitle((CFStep) i));
185 }
186
187 //____________________________________________________________________
188 AliUEHist::~AliUEHist()
189 {
190   // Destructor
191   
192   for (Int_t i=0; i<fkRegions; i++)
193   {
194     if (fTrackHist[i])
195     {
196       delete fTrackHist[i];
197       fTrackHist[i] = 0;
198     }
199   }
200      
201   if (fEventHist)
202   {
203     delete fEventHist;
204     fEventHist = 0;
205   }
206   
207   if (fCache)
208   {
209     delete fCache;
210     fCache = 0;
211   }
212 }
213
214 //____________________________________________________________________
215 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
216 {
217   // assigment operator
218
219   if (this != &c)
220     ((AliUEHist &) c).Copy(*this);
221
222   return *this;
223 }
224
225 //____________________________________________________________________
226 void AliUEHist::Copy(TObject& c) const
227 {
228   // copy function
229
230   AliUEHist& target = (AliUEHist &) c;
231
232   for (Int_t i=0; i<fkRegions; i++)
233     if (fTrackHist[i])
234       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
235
236   if (fEventHist)
237     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
238 }
239
240 //____________________________________________________________________
241 Long64_t AliUEHist::Merge(TCollection* list)
242 {
243   // Merge a list of AliUEHist objects with this (needed for
244   // PROOF). 
245   // Returns the number of merged objects (including this).
246
247   if (!list)
248     return 0;
249   
250   if (list->IsEmpty())
251     return 1;
252
253   TIterator* iter = list->MakeIterator();
254   TObject* obj;
255
256   // collections of objects
257   const Int_t kMaxLists = fkRegions+1;
258   TList** lists = new TList*[kMaxLists];
259   
260   for (Int_t i=0; i<kMaxLists; i++)
261     lists[i] = new TList;
262   
263   Int_t count = 0;
264   while ((obj = iter->Next())) {
265     
266     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
267     if (entry == 0) 
268       continue;
269
270     for (Int_t i=0; i<fkRegions; i++)
271       if (entry->fTrackHist[i])
272         lists[i]->Add(entry->fTrackHist[i]);
273     
274     lists[fkRegions]->Add(entry->fEventHist);
275
276     count++;
277   }
278   for (Int_t i=0; i<fkRegions; i++)
279     if (fTrackHist[i])
280       fTrackHist[i]->Merge(lists[i]);
281   
282   fEventHist->Merge(lists[fkRegions]);
283
284   for (Int_t i=0; i<kMaxLists; i++)
285     delete lists[i];
286     
287   delete[] lists;
288
289   return count+1;
290 }
291
292 //____________________________________________________________________
293 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
294 {
295   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
296   
297   if (fEtaMax > fEtaMin)
298     grid->SetRangeUser(0, fEtaMin, fEtaMax);
299   if (fPtMax > fPtMin)
300     grid->SetRangeUser(1, fPtMin, fPtMax);
301 }  
302
303 //____________________________________________________________________
304 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
305 {
306   // resets all bin limits 
307   
308   for (Int_t i=0; i<grid->GetNVar(); i++)
309     if (grid->GetGrid()->GetAxis(i)->TestBit(TH1::kIsZoomed))
310       grid->SetRangeUser(i, 0, -1);
311 }
312   
313 //____________________________________________________________________
314 TH1D* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax)
315 {
316   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
317   //
318   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
319   // Histogram has to be deleted by the caller of the function
320   
321   // unzoom all axes
322   ResetBinLimits(fTrackHist[region]->GetGrid(step));
323   ResetBinLimits(fEventHist->GetGrid(step));
324   
325   SetBinLimits(fTrackHist[region]->GetGrid(step));
326     
327   TH1D* tracks = 0;
328   
329   if (ptLeadMin < 0)
330   {
331     tracks = fTrackHist[region]->ShowProjection(2, step);
332     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
333     if (fCombineMinMax && region == kMin)
334     {
335       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
336       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
337       
338       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
339       tracks->Add(tracks2);
340       
341       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
342     }
343       
344     // normalize to get a density (deta dphi)
345     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
346     Float_t phiRegion = TMath::TwoPi() / 3;
347     if (!fCombineMinMax && region == kMin)
348       phiRegion /= 2;
349     Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
350     //Printf("Normalization: %f", normalization);
351     tracks->Scale(1.0 / normalization);
352     
353     TH1D* events = fEventHist->ShowProjection(0, step);
354     tracks->Divide(events);
355   }
356   else
357   {
358     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, ptLeadMin, ptLeadMax);
359     tracks = fTrackHist[region]->GetGrid(step)->Project(4);
360     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
361     
362     // normalize to get a density (deta dphi)
363     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
364     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1) * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
365     //Printf("Normalization: %f", normalization);
366     tracks->Scale(1.0 / normalization);
367     
368     TH1D* events = fEventHist->ShowProjection(0, step);
369     Int_t nEvents = events->Integral(events->FindBin(ptLeadMin), events->FindBin(ptLeadMax));
370     if (nEvents > 0)
371       tracks->Scale(1.0 / nEvents);
372   }
373
374   ResetBinLimits(fTrackHist[region]->GetGrid(step));
375
376   return tracks;
377 }
378
379 //____________________________________________________________________
380 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
381 {
382   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
383   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
384   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
385   //
386   // if trackCorrection is 0, just copies content from step1 to step2
387   
388   for (Int_t region=0; region<fkRegions; region++)
389   {
390     if (!fTrackHist[region])
391       continue;
392       
393     THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
394     THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
395     
396     if (trackCorrection != 0)
397     {
398       if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
399         AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
400         
401       if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
402         AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
403     }
404     
405     // optimized implementation
406     for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
407     {
408       Int_t bins[5];
409       Double_t value = grid->GetBinContent(binIdx, bins);
410       Double_t error = grid->GetBinError(binIdx);
411       
412       if (trackCorrection != 0)
413       {
414         if (var2 < 0)
415         {
416           value *= trackCorrection->GetBinContent(bins[var1]);
417           error *= trackCorrection->GetBinContent(bins[var1]);
418         }
419         else
420         {
421           value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
422           error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
423         }
424       }
425       
426       target->SetBinContent(bins, value);
427       target->SetBinError(bins, error);
428     }
429   }
430 }
431
432 //____________________________________________________________________
433 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection, Int_t var)
434 {
435   // corrects from step1 to step2 by multiplying the events with eventCorrection
436   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
437   //
438   // if eventCorrection is 0, just copies content from step1 to step2
439   
440   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
441   AliCFGridSparse* target = fEventHist->GetGrid(step2);
442   
443   if (eventCorrection != 0 && grid->GetNBins(var) != eventCorrection->GetNbinsX())
444     AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var), eventCorrection->GetNbinsX()));
445   
446   Int_t bins[2];
447   for (Int_t x = 1; x <= grid->GetNBins(0); x++)
448   {
449     bins[0] = x;
450     for (Int_t y = 1; y <= grid->GetNBins(1); y++)
451     {
452       bins[1] = y;
453       
454       Double_t value = grid->GetElement(bins);
455       if (value != 0)
456       {
457         Double_t error = grid->GetElementError(bins);
458         
459         if (eventCorrection != 0)
460         {
461           value *= eventCorrection->GetBinContent(bins[var]);
462           error *= eventCorrection->GetBinContent(bins[var]);
463         }
464         
465         target->SetElement(bins, value);
466         target->SetElementError(bins, error);
467       }
468     }
469   }
470 }
471
472 //____________________________________________________________________
473 void AliUEHist::Correct(AliUEHist* corrections)
474 {
475   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
476   //
477   // in this object the data is expected in the step kCFStepReconstructed
478   
479   // ---- track level
480   
481   // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
482   // extracted as function of leading pT
483   //TH1D* leadingBias = (TH1D*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, "z"); // from MC
484   TH1D* leadingBias = (TH1D*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, "z");              // from data
485   CorrectTracks(kCFStepReconstructed, kCFStepTracked, leadingBias, 2);
486   CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, 0);
487   delete leadingBias;
488   
489   // correct with kCFStepTracked --> kCFStepTrackedOnlyPrim
490   TH2D* contamination = corrections->GetTrackingContamination();
491   CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
492   CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
493   delete contamination;
494   
495   // correct with kCFStepTracked --> kCFStepAnaTopology
496   TH2D* trackingCorrection = corrections->GetTrackingCorrection();
497   CorrectTracks(kCFStepTracked, kCFStepAnaTopology, trackingCorrection, 0, 1);
498   CorrectEvents(kCFStepTracked, kCFStepAnaTopology, 0, 0);
499   delete trackingCorrection;
500   
501   // copy 
502   CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
503   CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
504   
505   // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
506   // practically independent of low pT cut 
507   TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
508   
509   // convert stage from true multiplicity to observed multiplicity by simple conversion factor
510   TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
511   vertexCorrectionObs->Reset();
512   
513   for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
514     vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(1.0 / 0.7 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i)));
515  
516 /*  new TCanvas;
517   vertexCorrection->DrawCopy();
518   vertexCorrectionObs->SetLineColor(2);
519   vertexCorrectionObs->DrawCopy("same");*/
520   
521   CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
522   CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
523   delete vertexCorrectionObs;
524   delete vertexCorrection;
525   
526   // copy 
527   CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
528   CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
529 }
530
531 //____________________________________________________________________
532 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2)
533 {
534   // creates a track-level efficiency by dividing step2 by step1
535   // projected to axis1 and axis2 (optional if >= 0)
536   
537   // integrate over regions
538   // cache it for efficiency (usually more than one efficiency is requested)
539   if (!fCache)
540   {
541     fCache = (AliCFContainer*) fTrackHist[0]->Clone();
542     for (Int_t i = 1; i < fkRegions; i++)
543       if (fTrackHist[i])
544         fCache->Add(fTrackHist[i]);
545   }
546       
547   // reset all limits and set the right ones except those in axis1 and axis2
548   ResetBinLimits(fCache->GetGrid(step1));
549   ResetBinLimits(fCache->GetGrid(step2));
550   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0)
551   {
552     fCache->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
553     fCache->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
554   }
555   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1)
556   {
557     fCache->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
558     fCache->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
559   }
560   
561   TH1* measured = 0;
562   TH1* generated = 0;
563     
564   if (axis2 > 0)
565   {
566     generated = fCache->Project(axis1, axis2, step1);
567     measured = fCache->Project(axis1, axis2, step2);
568   }
569   else
570   {
571     generated = fCache->Project(axis1, step1);
572     measured = fCache->Project(axis1, step2);
573   }
574   
575   measured->Divide(generated);
576   
577   delete generated;
578   
579   return measured;
580 }
581
582 //____________________________________________________________________
583 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
584 {
585   // creates a event-level efficiency by dividing step2 by step1
586   // projected to axis1 and axis2 (optional if >= 0)
587   
588   if (ptLeadMax > ptLeadMin)
589   {
590     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
591     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
592   }
593   
594   TH1* measured = 0;
595   TH1* generated = 0;
596     
597   if (axis2 > 0)
598   {
599     generated = fEventHist->Project(axis1, axis2, step1);
600     measured = fEventHist->Project(axis1, axis2, step2);
601   }
602   else
603   {
604     generated = fEventHist->Project(axis1, step1);
605     measured = fEventHist->Project(axis1, step2);
606   }
607   
608   measured->Divide(generated);
609   
610   delete generated;
611   
612   if (ptLeadMax > ptLeadMin)
613   {
614     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
615     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
616   }
617   
618   return measured;
619 }
620
621 //____________________________________________________________________
622 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
623 {
624   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2 
625   // where the matching is done of the z axis of hist1 with the x axis of hist2
626   
627   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
628     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
629   
630   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
631   {
632     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
633     {
634       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
635       {
636         if (hist2->GetBinContent(z) > 0)
637         {
638           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
639           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
640         }
641         else
642         {
643           hist1->SetBinContent(x, y, z, 0);
644           hist1->SetBinError(x, y, z, 0);
645         }
646       }
647     }
648   }
649 }  
650
651 //____________________________________________________________________
652 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, const char* axis, Float_t leadPtMin, Float_t leadPtMax)
653 {
654   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
655   // done by weighting the track-level distribution with the number of events as function of leading pT
656   // and then calculating the ratio between the distributions
657   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
658   //   no projection is done if axis == 0
659   
660   // integrate over regions
661   AliCFContainer* tmp = (AliCFContainer*) fTrackHist[0]->Clone();
662   for (Int_t i = 1; i < fkRegions; i++)
663     if (fTrackHist[i])
664       tmp->Add(fTrackHist[i]);
665   
666   TH1D* events1 = fEventHist->Project(0, step1);
667   TH3D* hist1 = tmp->Project(0, 1, 2, step1);
668   WeightHistogram(hist1, events1);
669   
670   TH1D* events2 = fEventHist->Project(0, step2);
671   TH3D* hist2 = tmp->Project(0, 1, 2, step2);
672   WeightHistogram(hist2, events2);
673   
674   TH1* generated = hist1;
675   TH1* measured = hist2;
676   
677   if (axis)
678   {
679     if (leadPtMax > leadPtMin)
680     {
681       hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
682       hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
683     }
684     
685     if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
686     {
687       hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
688       hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
689     }
690     if (fPtMax > fPtMin && !TString(axis).Contains("y"))
691     {
692       hist1->GetYaxis()->SetRangeUser(fPtMin, fPtMax);
693       hist2->GetYaxis()->SetRangeUser(fPtMin, fPtMax);
694     }
695     
696     generated = hist1->Project3D(axis);
697     measured  = hist2->Project3D(axis);
698     
699     // delete hists here if projection has been done
700     delete hist1;
701     delete hist2;
702   }
703   
704   measured->Divide(generated);
705   
706   delete events1;
707   delete events2;
708   delete generated;
709   delete tmp;
710   
711   return measured;
712 }
713
714 //____________________________________________________________________
715 TH2D* AliUEHist::GetTrackingEfficiency()
716 {
717   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
718   // integrates over the regions and all other variables than pT and eta to increase the statistics
719   //
720   // returned histogram has to be deleted by the user
721
722   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
723 }
724   
725 //____________________________________________________________________
726 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
727 {
728   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
729   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
730
731   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
732 }
733
734 //____________________________________________________________________
735 TH2D* AliUEHist::GetTrackingCorrection()
736 {
737   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
738   // integrates over the regions and all other variables than pT and eta to increase the statistics
739   //
740   // returned histogram has to be deleted by the user
741
742   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
743 }
744   
745 //____________________________________________________________________
746 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
747 {
748   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
749   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
750
751   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
752 }
753
754 //____________________________________________________________________
755 TH2D* AliUEHist::GetTrackingContamination()
756 {
757   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
758   // integrates over the regions and all other variables than pT and eta to increase the statistics
759   //
760   // returned histogram has to be deleted by the user
761
762   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
763 }
764   
765 //____________________________________________________________________
766 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
767 {
768   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
769   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
770
771   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
772 }
773
774 //____________________________________________________________________
775 const char* AliUEHist::GetRegionTitle(Region region)
776 {
777   // returns the name of the given region
778   
779   switch (region)
780   {
781     case kToward:
782       return "Towards";
783     case kAway:
784       return "Away";
785     case kMin:
786       return (fCombineMinMax) ? "Transverse" : "Min";
787     case kMax:
788       return "Max";
789   }
790   
791   return 0;
792 }
793   
794 //____________________________________________________________________
795 const char* AliUEHist::GetStepTitle(CFStep step)
796 {
797   // returns the name of the given step
798   
799   switch (step)
800   {
801     case kCFStepAll:
802       return "All events";
803     case kCFStepTriggered:
804       return "Triggered";
805     case kCFStepVertex:
806       return "Primary Vertex";
807     case kCFStepAnaTopology:
808       return "Required analysis topology";
809     case kCFStepTrackedOnlyPrim:
810       return "Tracked (matched MC, only primaries)";
811     case kCFStepTracked:
812       return "Tracked (matched MC, all)";
813     case kCFStepReconstructed:
814       return "Reconstructed";
815     case kCFStepRealLeading:
816       return "Correct leading particle identified";
817     case kCFStepBiasStudy:
818       return "Bias study applying tracking efficiency";
819     case kCFStepBiasStudy2:
820       return "Bias study applying tracking efficiency in two steps";
821   }
822   
823   return 0;
824 }