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