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