moving filling of V0 histogram later
[u/mrichter/AliRoot.git] / PWGCF / Correlations / Base / 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 #include "AliTHn.h"
38 #include "THn.h"
39
40 ClassImp(AliUEHist)
41
42 const Int_t AliUEHist::fgkCFSteps = 11;
43
44 AliUEHist::AliUEHist(const char* reqHist, const char* binning) : 
45   TObject(),
46   fkRegions(4),
47   fEventHist(0),
48   fTrackHistEfficiency(0),
49   fFakePt(0),
50   fEtaMin(0),
51   fEtaMax(0),
52   fPtMin(0),
53   fPtMax(0),
54   fPartSpecies(-1),
55   fCentralityMin(0),
56   fCentralityMax(0),
57   fZVtxMin(0),
58   fZVtxMax(0),
59   fContaminationEnhancement(0),
60   fCombineMinMax(0),
61   fTrackEtaCut(0),
62   fWeightPerEvent(0),
63   fCache(0),
64   fGetMultCacheOn(kFALSE),
65   fGetMultCache(0),
66   fHistogramType(reqHist)
67 {
68   // Constructor
69     
70   for (UInt_t i=0; i<fkRegions; i++)
71     fTrackHist[i] = 0;
72     
73   if (strlen(reqHist) == 0)
74     return;
75   
76   Printf("Creating AliUEHist with %s (binning: %s)", reqHist, binning);
77     
78   AliLog::SetClassDebugLevel("AliCFContainer", -1);
79   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
80     
81   const char* title = "";
82     
83   // track level
84   Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
85   Int_t iTrackBin[6];
86   Double_t* trackBins[6];
87   const char* trackAxisTitle[6];
88   
89   // eta
90   Int_t nEtaBins = -1;
91   Double_t* etaBins = GetBinning(binning, "eta", nEtaBins);
92   const char* etaTitle = "#eta";
93   
94   iTrackBin[0] = nEtaBins;
95   trackBins[0] = etaBins;
96   trackAxisTitle[0] = etaTitle;
97   
98   // delta eta
99   Int_t nDeltaEtaBins = -1;
100   Double_t* deltaEtaBins = GetBinning(binning, "delta_eta", nDeltaEtaBins);
101
102   // pT
103   trackBins[1] = GetBinning(binning, "p_t_assoc", iTrackBin[1]);
104   trackAxisTitle[1] = "p_{T} (GeV/c)";
105   
106   // pT, fine
107   Int_t npTBinsFine = -1;
108   Double_t* pTBinsFine = GetBinning(binning, "p_t_eff", npTBinsFine);
109
110   // pT,lead binning 1
111   Int_t nLeadingpTBins = -1;
112   Double_t* leadingpTBins = GetBinning(binning, "p_t_leading", nLeadingpTBins);
113
114   // pT,lead binning 2
115   Int_t nLeadingpTBins2 = -1;
116   Double_t* leadingpTBins2 = GetBinning(binning, "p_t_leading_course", nLeadingpTBins2);
117   
118   // phi,lead
119   Int_t nLeadingPhiBins = -1;
120   Double_t* leadingPhiBins = GetBinning(binning, "delta_phi", nLeadingPhiBins);
121   
122   trackBins[3] = GetBinning(binning, "multiplicity", iTrackBin[3]);
123   trackAxisTitle[3] = "multiplicity";
124   
125   // particle species
126   const Int_t kNSpeciesBins = 4; // pi, K, p, rest
127   Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
128   
129   // vtx-z axis
130   const char* vertexTitle = "z-vtx (cm)";
131   Int_t nVertexBins = -1;
132   Double_t* vertexBins = GetBinning(binning, "vertex", nVertexBins);
133   Int_t nVertexBinsEff = -1;
134   Double_t* vertexBinsEff = GetBinning(binning, "vertex_eff", nVertexBinsEff);
135   
136   Int_t useVtxAxis = 0;
137   
138   // selection depending on requested histogram
139   Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
140   if (strcmp(reqHist, "NumberDensitypT") == 0)
141   {
142     axis = 0;
143     title = "d^{2}N_{ch}/d#varphid#eta";
144   }
145   else if (strcmp(reqHist, "NumberDensityPhi") == 0)
146   {
147     axis = 1;
148     title = "d^{2}N_{ch}/d#varphid#eta";
149   }
150   else if (TString(reqHist).BeginsWith("NumberDensityPhiCentrality"))
151   {
152     if (TString(reqHist).Contains("Vtx"))
153       useVtxAxis = 1;
154     
155     reqHist = "NumberDensityPhiCentrality";
156     fHistogramType = reqHist;
157     axis = 2;
158     title = "d^{2}N_{ch}/d#varphid#eta";
159   }
160   else if (strcmp(reqHist, "SumpT") == 0)
161   {
162     axis = 0;
163     title = "d^{2}#Sigma p_{T}/d#varphid#eta";
164   }
165   else
166     AliFatal(Form("Invalid histogram requested: %s", reqHist));
167   
168   UInt_t initRegions = fkRegions;
169   
170   if (axis == 0)
171   {
172     trackBins[2] = leadingpTBins;
173     iTrackBin[2] = nLeadingpTBins;
174     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
175     
176   }
177   else if (axis == 1)
178   {
179     nTrackVars = 5;
180     initRegions = 1;
181   
182     iTrackBin[2] = nLeadingpTBins2;
183     trackBins[2] = leadingpTBins2;
184     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
185     
186     iTrackBin[4] = nLeadingPhiBins;
187     trackBins[4] = leadingPhiBins;
188     trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
189   }
190   else if (axis == 2)
191   {
192     nTrackVars = 5;
193     initRegions = 1;
194   
195     iTrackBin[0] = nDeltaEtaBins;
196     trackBins[0] = deltaEtaBins;
197     trackAxisTitle[0] = "#Delta#eta";
198   
199     iTrackBin[2] = nLeadingpTBins2;
200     trackBins[2] = leadingpTBins2;
201     trackAxisTitle[2] = "leading p_{T} (GeV/c)";
202     
203     trackAxisTitle[3] = "centrality";
204   
205     iTrackBin[4] = nLeadingPhiBins;
206     trackBins[4] = leadingPhiBins;
207     trackAxisTitle[4] = "#Delta#varphi (rad)";
208
209     if (useVtxAxis > 0)
210     {
211       nTrackVars = 6;
212       iTrackBin[5] = nVertexBins;
213       trackBins[5] = vertexBins;
214       trackAxisTitle[5] = vertexTitle;
215     }
216   }
217     
218   for (UInt_t i=0; i<initRegions; i++)
219   {
220     if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
221       fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
222     else
223       fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
224     
225     for (Int_t j=0; j<nTrackVars; j++)
226     {
227       fTrackHist[i]->SetBinLimits(j, trackBins[j]);
228       fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
229     }
230     
231     SetStepNames(fTrackHist[i]);
232   }
233   
234   // event level
235   Int_t nEventVars = 2;
236   Int_t iEventBin[3];
237
238   // track 3rd and 4th axis --> event 1st and 2nd axis
239   iEventBin[0] = iTrackBin[2];
240   iEventBin[1] = iTrackBin[3];
241   
242   // plus track 5th axis (in certain cases)
243   if (axis == 2 && useVtxAxis)
244   {
245     nEventVars = 3;
246     iEventBin[2] = iTrackBin[5];
247   }
248   
249   fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
250   
251   fEventHist->SetBinLimits(0, trackBins[2]);
252   fEventHist->SetVarTitle(0, trackAxisTitle[2]);
253   
254   fEventHist->SetBinLimits(1, trackBins[3]);
255   fEventHist->SetVarTitle(1, trackAxisTitle[3]);
256   
257   if (axis == 2 && useVtxAxis)
258   {
259     fEventHist->SetBinLimits(2, trackBins[5]);
260     fEventHist->SetVarTitle(2, trackAxisTitle[5]);
261   }
262
263   SetStepNames(fEventHist);
264   
265   iTrackBin[0] = nEtaBins;
266   iTrackBin[1] = npTBinsFine;
267   iTrackBin[2] = kNSpeciesBins;
268   iTrackBin[4] = nVertexBinsEff;
269
270   fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 6, 5, iTrackBin);
271   fTrackHistEfficiency->SetBinLimits(0, etaBins);
272   fTrackHistEfficiency->SetVarTitle(0, etaTitle);
273   fTrackHistEfficiency->SetBinLimits(1, pTBinsFine);
274   fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
275   fTrackHistEfficiency->SetBinLimits(2, speciesBins);
276   fTrackHistEfficiency->SetVarTitle(2, "particle species");
277   fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
278   fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
279   fTrackHistEfficiency->SetBinLimits(4, vertexBinsEff);
280   fTrackHistEfficiency->SetVarTitle(4, vertexTitle);
281
282   fFakePt = new TH3F("fFakePt","fFakePt;p_{T,rec};p_{T};centrality", 200, 0, 20, 200, 0, 20, 20, 0, 100);
283   
284   delete[] deltaEtaBins;
285   delete[] pTBinsFine;
286   delete[] leadingpTBins;
287   delete[] leadingpTBins2;
288   delete[] leadingPhiBins;
289   delete[] vertexBins;
290   delete[] vertexBinsEff;
291 }
292
293 Double_t* AliUEHist::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
294 {
295   // takes the binning from <configuration> identified by <tag>
296   // configuration syntax example:
297   // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
298   // phi: .....
299   //
300   // returns bin edges which have to be deleted by the caller
301   
302   TString config(configuration);
303   TObjArray* lines = config.Tokenize("\n");
304   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
305   {
306     TString line(lines->At(i)->GetName());
307     if (line.BeginsWith(TString(tag) + ":"))
308     {
309       line.Remove(0, strlen(tag) + 1);
310       line.ReplaceAll(" ", "");
311       TObjArray* binning = line.Tokenize(",");
312       Double_t* bins = new Double_t[binning->GetEntriesFast()];
313       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
314         bins[j] = TString(binning->At(j)->GetName()).Atof();
315       
316       nBins = binning->GetEntriesFast() - 1;
317
318       delete binning;
319       delete lines;
320       return bins;
321     }
322   }
323   
324   delete lines;
325   AliFatal(Form("Tag %s not found in %s", tag, configuration));
326   return 0;
327 }
328
329 //_____________________________________________________________________________
330 AliUEHist::AliUEHist(const AliUEHist &c) :
331   TObject(),
332   fkRegions(4),
333   fEventHist(0),
334   fTrackHistEfficiency(0),
335   fFakePt(0),
336   fEtaMin(0),
337   fEtaMax(0),
338   fPtMin(0),
339   fPtMax(0),
340   fPartSpecies(-1),
341   fCentralityMin(0),
342   fCentralityMax(0),
343   fZVtxMin(0),
344   fZVtxMax(0),
345   fContaminationEnhancement(0),
346   fCombineMinMax(0),
347   fTrackEtaCut(0),
348   fWeightPerEvent(0),
349   fCache(0),
350   fGetMultCacheOn(kFALSE),
351   fGetMultCache(0),
352   fHistogramType()
353 {
354   //
355   // AliUEHist copy constructor
356   //
357
358   ((AliUEHist &) c).Copy(*this);
359 }
360
361 //____________________________________________________________________
362 void AliUEHist::SetStepNames(AliCFContainer* container)
363 {
364   // sets the names of the correction steps
365   
366   for (Int_t i=0; i<fgkCFSteps; i++)
367     container->SetStepTitle(i, GetStepTitle((CFStep) i));
368 }
369
370 //____________________________________________________________________
371 AliUEHist::~AliUEHist()
372 {
373   // Destructor
374   
375   for (UInt_t i=0; i<fkRegions; i++)
376   {
377     if (fTrackHist[i])
378     {
379       delete fTrackHist[i];
380       fTrackHist[i] = 0;
381     }
382   }
383      
384   if (fEventHist)
385   {
386     delete fEventHist;
387     fEventHist = 0;
388   }
389   
390   if (fTrackHistEfficiency)
391   {
392     delete fTrackHistEfficiency;
393     fTrackHistEfficiency = 0;
394   }
395
396   if (fFakePt)
397   {
398     delete fFakePt;
399     fFakePt = 0;
400   }
401
402   if (fCache)
403   {
404     delete fCache;
405     fCache = 0;
406   }
407 }
408
409 //____________________________________________________________________
410 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
411 {
412   // assigment operator
413
414   if (this != &c)
415     ((AliUEHist &) c).Copy(*this);
416
417   return *this;
418 }
419
420 //____________________________________________________________________
421 void AliUEHist::Copy(TObject& c) const
422 {
423   // copy function
424
425   AliUEHist& target = (AliUEHist &) c;
426
427   for (UInt_t i=0; i<fkRegions; i++)
428     if (fTrackHist[i])
429       target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
430
431   if (fEventHist)
432     target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
433   
434   if (fTrackHistEfficiency)
435     target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
436     
437   if (fFakePt)
438     target.fFakePt = dynamic_cast<TH3F*> (fFakePt->Clone());
439
440   target.fEtaMin = fEtaMin;
441   target.fEtaMax = fEtaMax;
442   target.fPtMin = fPtMin;
443   target.fPtMax = fPtMax;
444   target.fPartSpecies = fPartSpecies;
445   target.fCentralityMin = fCentralityMin;
446   target.fCentralityMax = fCentralityMax;
447   target.fZVtxMin = fZVtxMin;
448   target.fZVtxMax = fZVtxMax;
449   
450   if (fContaminationEnhancement)
451     target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
452     
453   target.fCombineMinMax = fCombineMinMax;
454   target.fTrackEtaCut = fTrackEtaCut;
455   target.fWeightPerEvent = fWeightPerEvent;
456   target.fHistogramType = fHistogramType;
457 }
458
459 //____________________________________________________________________
460 Long64_t AliUEHist::Merge(TCollection* list)
461 {
462   // Merge a list of AliUEHist objects with this (needed for
463   // PROOF). 
464   // Returns the number of merged objects (including this).
465
466   if (!list)
467     return 0;
468   
469   if (list->IsEmpty())
470     return 1;
471
472   TIterator* iter = list->MakeIterator();
473   TObject* obj;
474
475   // collections of objects
476   const UInt_t kMaxLists = fkRegions+3;
477   TList** lists = new TList*[kMaxLists];
478   
479   for (UInt_t i=0; i<kMaxLists; i++)
480     lists[i] = new TList;
481   
482   Int_t count = 0;
483   while ((obj = iter->Next())) {
484     
485     AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
486     if (entry == 0) 
487       continue;
488
489     for (UInt_t i=0; i<fkRegions; i++)
490       if (entry->fTrackHist[i])
491         lists[i]->Add(entry->fTrackHist[i]);
492     
493     lists[fkRegions]->Add(entry->fEventHist);
494     lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
495     if (entry->fFakePt)
496       lists[fkRegions+2]->Add(entry->fFakePt);
497
498     count++;
499   }
500   for (UInt_t i=0; i<fkRegions; i++)
501     if (fTrackHist[i])
502       fTrackHist[i]->Merge(lists[i]);
503   
504   fEventHist->Merge(lists[fkRegions]);
505   fTrackHistEfficiency->Merge(lists[fkRegions+1]);
506   if (fFakePt)
507     fFakePt->Merge(lists[fkRegions+2]);
508
509   for (UInt_t i=0; i<kMaxLists; i++)
510     delete lists[i];
511     
512   delete[] lists;
513
514   return count+1;
515 }
516
517 //____________________________________________________________________
518 void AliUEHist::SetBinLimits(THnBase* grid)
519 {
520   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
521   
522   if (fEtaMax > fEtaMin)
523     grid->GetAxis(0)->SetRangeUser(fEtaMin, fEtaMax);
524   if (fPtMax > fPtMin)
525     grid->GetAxis(1)->SetRangeUser(fPtMin, fPtMax);
526 }  
527
528 //____________________________________________________________________
529 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
530 {
531   // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
532   
533   SetBinLimits(grid->GetGrid());
534 }
535
536 //____________________________________________________________________
537 void AliUEHist::ResetBinLimits(THnBase* grid)
538 {
539   // resets all bin limits 
540
541   for (Int_t i=0; i<grid->GetNdimensions(); i++)
542     if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange))
543       grid->GetAxis(i)->SetRangeUser(0, -1);
544 }
545
546 //____________________________________________________________________
547 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
548 {
549   // resets all bin limits 
550   
551   ResetBinLimits(grid->GetGrid());
552 }
553
554 //____________________________________________________________________
555 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
556 {
557   // prints the number of empty bins in the track end event histograms in the given step
558   
559   Int_t binBegin[4];
560   Int_t binEnd[4];
561   
562   for (Int_t i=0; i<4; i++)
563   {
564     binBegin[i] = 1;
565     binEnd[i]   = fTrackHist[0]->GetNBins(i);
566   }
567   
568   if (fEtaMax > fEtaMin)
569   {
570     binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
571     binEnd[0]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
572   }
573   
574   if (fPtMax > fPtMin)
575   {
576     binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
577     binEnd[1]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
578   }
579   
580   if (ptLeadMax > ptLeadMin)
581   {
582     binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
583     binEnd[2]   = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
584   }
585   
586   // start from multiplicity 1
587   binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
588   
589   for (UInt_t region=0; region<fkRegions; region++)
590   {
591     Int_t total = 0;
592     Int_t count = 0;
593     Int_t vars[4];
594     
595     for (Int_t i=0; i<4; i++)
596       vars[i] = binBegin[i];
597       
598     AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
599     while (1)
600     {
601       if (grid->GetElement(vars) == 0)
602       {
603         Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f", 
604           grid->GetBinCenter(0, vars[0]), 
605           grid->GetBinCenter(1, vars[1]), 
606           grid->GetBinCenter(2, vars[2]), 
607           grid->GetBinCenter(3, vars[3])
608         );
609         count++;
610       }
611       
612       vars[3]++;
613       for (Int_t i=3; i>0; i--)
614       {
615         if (vars[i] == binEnd[i]+1)
616         {
617           vars[i] = binBegin[i];
618           vars[i-1]++;
619         }
620       }
621       
622       if (vars[0] == binEnd[0]+1)
623         break;
624       total++;
625     }
626   
627     Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
628   }
629 }  
630
631 //____________________________________________________________________
632 TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm, Long64_t* normEvents)
633 {
634   // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
635   //
636   // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
637   // Histogram has to be deleted by the caller of the function
638   //
639   // twoD: 0: 1D histogram as function of phi
640   //       1: 2D histogram, deltaphi, deltaeta
641   //       10: 1D histogram, within |deltaeta| < 0.8
642   //       11: 1D histogram, within 0.8 < |deltaeta| < 1.6
643   //
644   // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
645   //
646   // normEvents: if non-0 the number of events/trigger particles for the normalization is filled
647   
648   // unzoom all axes
649   ResetBinLimits(fTrackHist[region]->GetGrid(step));
650   ResetBinLimits(fEventHist->GetGrid(step));
651   
652   SetBinLimits(fTrackHist[region]->GetGrid(step));
653   
654   // z vtx
655   if (fZVtxMax > fZVtxMin)
656   {
657     Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
658     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
659     fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
660   }
661     
662   TH1D* tracks = 0;
663   
664   if (ptLeadMin < 0)
665   {
666     tracks = fTrackHist[region]->ShowProjection(2, step);
667     tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
668     if (fCombineMinMax && region == kMin)
669     {
670       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
671       SetBinLimits(fTrackHist[kMax]->GetGrid(step));
672       
673       TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
674       tracks->Add(tracks2);
675       
676       ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
677     }
678       
679     // normalize to get a density (deta dphi)
680     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
681     Float_t phiRegion = TMath::TwoPi() / 3;
682     if (!fCombineMinMax && region == kMin)
683       phiRegion /= 2;
684     Float_t normalization = phiRegion;
685     if (etaNorm)
686       normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
687     //Printf("Normalization: %f", normalization);
688     tracks->Scale(1.0 / normalization);
689     
690     TH1D* events = fEventHist->ShowProjection(0, step);
691     tracks->Divide(events);
692   }
693   else
694   {
695     if (multBinEnd >= multBinBegin)
696     {
697       Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
698       fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
699       fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
700     }
701     
702     Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
703     Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
704     
705     Printf("Using leading pT range %d --> %d", firstBin, lastBin);
706     
707     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
708     
709     if (twoD == 0)
710       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
711     else
712       tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
713       
714     Printf("Calculated histogram --> %f tracks", tracks->Integral());
715     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
716     
717     if (twoD == 10 || twoD == 11)
718     {
719       Float_t etaLimit = 1.0;
720       if (twoD == 10)
721       {
722         tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
723         
724         // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
725         tracks->Scale(1. / 0.75);
726       }
727       else if (twoD == 11)
728       {
729         TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
730         TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
731         
732         tracksTmp1->Add(tracksTmp2);
733         
734         delete tracks;
735         tracks = tracksTmp1;
736         delete tracksTmp2;
737         
738         tracks->Scale(1. / 0.25);
739       }
740     }
741     
742     // normalize to get a density (deta dphi)
743     // TODO normalization may be off for 2d histogram
744     Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
745     
746     if (etaNorm)
747     {
748       TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
749       if (strcmp(axis->GetTitle(), "#eta") == 0)
750       {
751         Printf("Normalizing using eta-axis range");
752         normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
753       }
754       else
755       {
756         Printf("Normalizing assuming accepted range of |eta| < 0.8");
757         normalization *= 0.8 * 2;
758       }
759     }
760     
761     //Printf("Normalization: %f", normalization);
762     tracks->Scale(1.0 / normalization);
763     
764     // NOTE fEventHist contains the number of events for the underlying event analysis and the number of trigger particles for the azimuthal correlation analysis. In the latter case the naming is therefore somewhat misleading!
765     TH1D* events = fEventHist->ShowProjection(0, step);
766     Long64_t nEvents = (Long64_t) events->Integral(firstBin, lastBin);
767     Printf("Calculated histogram --> %lld events", nEvents);
768     if (normEvents)
769       *normEvents = nEvents;
770       
771     if (nEvents > 0)
772       tracks->Scale(1.0 / nEvents);
773     
774     delete events;
775   }
776
777   ResetBinLimits(fTrackHist[region]->GetGrid(step));
778   ResetBinLimits(fEventHist->GetGrid(step));
779
780   return tracks;
781 }
782
783 //____________________________________________________________________
784 void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist)
785 {
786   // Calculates a 3d histogram with deltaphi, deltaeta, zvtx on track level and a 1d histogram on event level (as fct of zvtx)
787   // Histogram has to be deleted by the caller of the function
788   
789   if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
790     AliFatal("Histogram without vertex axis provided");
791
792   // unzoom all axes
793   ResetBinLimits(fTrackHist[region]->GetGrid(step));
794   ResetBinLimits(fEventHist->GetGrid(step));
795   
796   SetBinLimits(fTrackHist[region]->GetGrid(step));
797   
798   if (multBinEnd >= multBinBegin)
799   {
800     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
801     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
802     fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
803   }
804     
805   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
806   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
807   Printf("Using leading pT range %d --> %d", firstBin, lastBin);
808   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
809   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
810     
811   *trackHist = (TH3*) fTrackHist[region]->GetGrid(step)->Project(4, 0, 5);
812   *eventHist = (TH1*) fEventHist->GetGrid(step)->Project(2);
813
814   ResetBinLimits(fTrackHist[region]->GetGrid(step));
815   ResetBinLimits(fEventHist->GetGrid(step));
816 }
817
818 //____________________________________________________________________
819 void AliUEHist::GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist)
820 {
821   // Calculates a 4d histogram with deltaphi, deltaeta, zvtx, multiplicity on track level and 
822   // a 2d histogram on event level (as fct of zvtx, multiplicity)
823   // Histograms has to be deleted by the caller of the function
824   
825   if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
826     AliFatal("Histogram without vertex axis provided");
827
828   THnBase* sparse = fTrackHist[region]->GetGrid(step)->GetGrid();
829   if (fGetMultCacheOn)
830   {
831     if (!fGetMultCache)
832     {
833       fGetMultCache = ChangeToThn(sparse);
834       // should work but causes SEGV in ProjectionND below 
835     }
836     sparse = fGetMultCache;
837   }
838   
839   // unzoom all axes
840   ResetBinLimits(sparse);
841   ResetBinLimits(fEventHist->GetGrid(step));
842   
843   SetBinLimits(sparse);
844   
845   Int_t firstBin = sparse->GetAxis(2)->FindBin(ptLeadMin);
846   Int_t lastBin = sparse->GetAxis(2)->FindBin(ptLeadMax);
847   Printf("Using leading pT range %d --> %d", firstBin, lastBin);
848   sparse->GetAxis(2)->SetRange(firstBin, lastBin);
849   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
850     
851   Int_t dimensions[] = { 4, 0, 5, 3 };
852   THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E");
853   *eventHist = (TH2*) fEventHist->GetGrid(step)->Project(2, 1);
854   
855   ResetBinLimits(sparse);
856   ResetBinLimits(fEventHist->GetGrid(step));
857
858   // convert to THn 
859   *trackHist = ChangeToThn(tmpTrackHist);
860   delete tmpTrackHist;
861 }
862
863 //____________________________________________________________________
864 TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t normalizePerTrigger)
865 {
866   // Calls GetUEHist(...) for *each* vertex bin and multiplicity bin and performs a sum of ratios:
867   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
868   // where N is the total number of events/trigger particles and the subscript is the vertex/multiplicity bin
869   // where mixed is normalized such that the information about the number of pairs in same is kept
870   //
871   // returns a 2D histogram: deltaphi, deltaeta
872   //
873   // Parameters:
874   //   mixed: AliUEHist containing mixed event corresponding to this object
875   //   <other parameters> : check documentation of AliUEHist::GetUEHist
876   //  normalizePerTrigger: divide through number of triggers
877   
878   // do not add this hists to the directory
879   Bool_t oldStatus = TH1::AddDirectoryStatus();
880   TH1::AddDirectory(kFALSE);
881
882   TH2* totalTracks = 0;
883   
884   THnBase* trackSameAll = 0;
885   THnBase* trackMixedAll = 0;
886   TH2* eventSameAll = 0;
887   TH2* eventMixedAll = 0;
888   
889   Int_t totalEvents = 0;
890   Int_t nCorrelationFunctions = 0;
891   
892   GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackSameAll, &eventSameAll);
893   mixed->GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackMixedAll, &eventMixedAll);
894   
895 //   Printf("%f %f %f %f", trackSameAll->GetEntries(), eventSameAll->GetEntries(), trackMixedAll->GetEntries(), eventMixedAll->GetEntries());
896   
897 //   TH1* normParameters = new TH1F("normParameters", "", 100, 0, 2);
898   
899 //   trackSameAll->Dump();
900
901   TAxis* multAxis = trackSameAll->GetAxis(3);
902
903   if (multBinEnd < multBinBegin)
904   {
905     multBinBegin = 1;
906     multBinEnd = multAxis->GetNbins();
907   }
908   
909   for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++)
910   {
911     trackSameAll->GetAxis(3)->SetRange(multBin, multBin);
912     trackMixedAll->GetAxis(3)->SetRange(multBin, multBin);
913
914     // get mixed normalization correction factor: is independent of vertex bin if scaled with number of triggers
915     trackMixedAll->GetAxis(2)->SetRange(0, -1);
916     TH2* tracksMixed = trackMixedAll->Projection(1, 0, "E");
917 //     Printf("%f", tracksMixed->Integral());
918     Float_t binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1);
919     
920     // get mixed event normalization by assuming full acceptance at deta at 0 (integrate over dphi), excluding (0, 0)
921     Double_t mixedNormError = 0;
922     Double_t mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetYaxis()->FindBin(-0.01)-1, tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
923     Double_t mixedNormError2 = 0;
924     Double_t mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetYaxis()->FindBin(0.01)+1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2);
925     
926     if (mixedNormError == 0 || mixedNormError2 == 0)
927     {
928       Printf("ERROR: Skipping multiplicity %d because mixed event is empty %f %f %f %f", multBin, mixedNorm, mixedNormError, mixedNorm2, mixedNormError2);
929       continue;
930     }
931     
932     Int_t nBinsMixedNorm = (tracksMixed->GetYaxis()->FindBin(-0.01) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
933     mixedNorm /= nBinsMixedNorm;
934     mixedNormError /= nBinsMixedNorm;
935
936     Int_t nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetYaxis()->FindBin(0.01) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
937     mixedNorm2 /= nBinsMixedNorm2;
938     mixedNormError2 /= nBinsMixedNorm2;
939
940     mixedNorm = mixedNorm / mixedNormError / mixedNormError + mixedNorm2 / mixedNormError2 / mixedNormError2;
941     mixedNormError = TMath::Sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2));
942     mixedNorm *= mixedNormError * mixedNormError;
943     
944 /*    Double_t mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
945     Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);*/
946
947     delete tracksMixed;
948     
949     Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin);
950 //     Printf("%f +- %f | %f | %f", mixedNorm, mixedNormError, triggers, mixedNorm / triggers);
951     if (triggers <= 0)
952     {
953       Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
954       continue;
955     }
956     
957     mixedNorm /= triggers;
958     mixedNormError /= triggers;
959     
960     if (mixedNorm <= 0)
961     {
962       Printf("ERROR: Skipping multiplicity %d because mixed event is empty at (0,0)", multBin);
963       continue;
964     }
965     
966     // finite bin correction
967     if (fTrackEtaCut > 0)
968     {
969       Double_t finiteBinCorrection = -1.0 / (2*fTrackEtaCut) * binWidthEta / 2 + 1;
970       Printf("Finite bin correction: %f", finiteBinCorrection);
971       mixedNorm /= finiteBinCorrection;
972       mixedNormError /= finiteBinCorrection;
973     }
974     else
975     {
976       Printf("ERROR: fTrackEtaCut not set. Finite bin correction cannot be applied. Continuing anyway...");
977     }
978     
979 //     Printf("Norm: %f +- %f", mixedNorm, mixedNormError);
980
981 //     normParameters->Fill(mixedNorm);
982       
983     TAxis* vertexAxis = trackSameAll->GetAxis(2);
984     for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
985 //     for (Int_t vertexBin = 5; vertexBin <= 6; vertexBin++)
986     {
987       trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
988       trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
989
990       TH2* tracksSame = trackSameAll->Projection(1, 0, "E");
991       tracksMixed = trackMixedAll->Projection(1, 0, "E");
992       
993       // asssume flat in dphi, gain in statistics
994       //     TH1* histMixedproj = mixedTwoD->ProjectionY();
995       //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
996       //     
997       //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
998       //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
999       //        mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1000
1001       //       delete histMixedproj;
1002       
1003       Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin);
1004       if (triggers2 <= 0)
1005       {
1006         Printf("ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin);
1007       }
1008       else
1009       {
1010         tracksMixed->Scale(1.0 / triggers2 / mixedNorm);
1011         // tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
1012           
1013 //      new TCanvas; tracksSame->DrawClone("SURF1");
1014 //      new TCanvas; tracksMixed->DrawClone("SURF1");
1015
1016         // some code to judge the relative contribution of the different correlation functions to the overall uncertainty
1017         Double_t sums[] = { 0, 0, 0 };
1018         Double_t errors[] = { 0, 0, 0 };
1019
1020         for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1021           for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1022           {
1023             sums[0] += tracksSame->GetBinContent(x, y);
1024             errors[0] += tracksSame->GetBinError(x, y);
1025             sums[1] += tracksMixed->GetBinContent(x, y);
1026             errors[1] += tracksMixed->GetBinError(x, y);
1027           }
1028
1029         tracksSame->Divide(tracksMixed);
1030           
1031         for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1032           for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1033           {
1034             sums[2] += tracksSame->GetBinContent(x, y);
1035             errors[2] += tracksSame->GetBinError(x, y);
1036           }
1037           
1038         for (Int_t x=0; x<3; x++)
1039           if (sums[x] > 0)
1040             errors[x] /= sums[x];
1041           
1042         Printf("The correlation function %d %d has uncertainties %f %f %f (Ratio S/M %f)", multBin, vertexBin, errors[0], errors[1], errors[2], (errors[1] > 0) ? errors[0] / errors[1] : -1);
1043         // code to draw contributions
1044         /*
1045         TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
1046         proj->SetTitle(Form("Bin %d", vertexBin));
1047         proj->SetLineColor(vertexBin);
1048         proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
1049         */
1050         
1051         if (!totalTracks)
1052           totalTracks = (TH2*) tracksSame->Clone("totalTracks");
1053         else
1054           totalTracks->Add(tracksSame);
1055
1056         totalEvents += eventSameAll->GetBinContent(vertexBin, multBin);
1057         
1058 //      new TCanvas; tracksMixed->DrawCopy("SURF1");
1059       }
1060
1061       delete tracksSame;
1062       delete tracksMixed;
1063       
1064       nCorrelationFunctions++;
1065     }
1066   }
1067
1068   if (totalTracks) {
1069     Double_t sums[] = { 0, 0, 0 };
1070     Double_t errors[] = { 0, 0, 0 };
1071
1072     for (Int_t x=1; x<=totalTracks->GetNbinsX(); x++)
1073       for (Int_t y=1; y<=totalTracks->GetNbinsY(); y++)
1074       {
1075         sums[0] += totalTracks->GetBinContent(x, y);
1076         errors[0] += totalTracks->GetBinError(x, y);
1077       }
1078     if (sums[0] > 0)
1079       errors[0] /= sums[0];
1080     
1081     if (normalizePerTrigger)
1082     {
1083       Printf("Dividing %f tracks by %d events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]);
1084       if (totalEvents > 0)
1085         totalTracks->Scale(1.0 / totalEvents);
1086     }
1087   
1088     // normalizate to dphi width
1089     Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1);
1090     totalTracks->Scale(1.0 / normalization);
1091   }
1092   
1093   delete trackSameAll;
1094   delete trackMixedAll;
1095   delete eventSameAll;
1096   delete eventMixedAll;
1097   
1098 //   new TCanvas; normParameters->Draw();
1099   
1100   TH1::AddDirectory(oldStatus);
1101   
1102   return totalTracks;
1103 }
1104
1105 TH1* AliUEHist::GetTriggersAsFunctionOfMultiplicity(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
1106 {
1107   // returns the distribution of triggers as function of centrality/multiplicity
1108
1109   ResetBinLimits(fEventHist->GetGrid(step));
1110   
1111   Int_t firstBin = fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(ptLeadMin);
1112   Int_t lastBin = fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(ptLeadMax);
1113   Printf("Using pT range %d --> %d", firstBin, lastBin);
1114   fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
1115   
1116   TH1* eventHist = fEventHist->GetGrid(step)->Project(1);
1117
1118   ResetBinLimits(fEventHist->GetGrid(step));
1119   
1120   return eventHist;  
1121 }
1122
1123 /*
1124 TH2* AliUEHist::GetPtDistInPhiRegion(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiBegin, Float_t phiEnd)
1125 {
1126   // Returns pT,assoc distribution in the given pt,trig, multiplicity, phi region
1127   // Does not use sum of ratios for mixed event correction (TODO to be improved)
1128   // returns a 2D histogram: deltaphi, deltaeta
1129   //
1130   // Parameters:
1131   //   mixed: AliUEHist containing mixed event corresponding to this object
1132 */
1133
1134 //____________________________________________________________________
1135 TH2* AliUEHist::GetSumOfRatios(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm, Bool_t useVertexBins)
1136 {
1137   // Calls GetUEHist(...) for *each* multiplicity bin and performs a sum of ratios:
1138   // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
1139   // where N is the total number of events/trigger particles and the subscript is the multiplicity bin
1140   //
1141   // Can only be used for the 2D histogram at present
1142   //
1143   // Parameters:
1144   //   mixed: AliUEHist containing mixed event corresponding to this object
1145   //   <other parameters> : check documentation of AliUEHist::GetUEHist
1146   
1147   TH2* totalTracks = 0;
1148   Int_t totalEvents = 0;
1149   
1150   Int_t vertexBin = 1;
1151   TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
1152   if (useVertexBins && !vertexAxis)
1153   {
1154     Printf("Vertex axis requested but not available");
1155     return 0;
1156   }
1157   
1158   // vertex bin loop
1159   while (1)
1160   {
1161     if (useVertexBins)
1162     {
1163       SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1164       mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1165       vertexBin++;
1166     }
1167     
1168     // multiplicity loop
1169     Int_t multIter = multBinBegin;
1170     while (1)
1171     {
1172       Int_t multBinBeginLocal = multBinBegin;
1173       Int_t multBinEndLocal = multBinEnd;
1174       
1175       if (multBinEnd >= multBinBegin)
1176       {
1177         multBinBeginLocal = multIter;
1178         multBinEndLocal = multIter;
1179         multIter++;
1180       }
1181         
1182       Long64_t nEvents = 0;
1183       TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
1184       // undo normalization
1185       tracks->Scale(nEvents);
1186       totalEvents += nEvents;
1187       
1188       TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
1189       
1190       // asssume flat in dphi, gain in statistics
1191   //     TH1* histMixedproj = mixedTwoD->ProjectionY();
1192   //     histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
1193   //     
1194   //     for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
1195   //       for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
1196   //    mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1197
1198   //       delete histMixedproj;
1199   
1200       // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
1201 /*      Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
1202       mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
1203       tracks->Scale(mixedNorm);*/
1204       
1205       tracks->Scale(mixedTwoD->Integral() / tracks->Integral());
1206
1207       tracks->Divide(mixedTwoD);
1208       
1209       delete mixedTwoD;
1210       
1211       if (!totalTracks)
1212         totalTracks = tracks;
1213       else
1214       {
1215         totalTracks->Add(tracks);
1216         delete tracks;
1217       }
1218
1219       if (multIter > multBinEnd)
1220         break;
1221     }
1222     
1223     if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
1224       break;
1225   }
1226
1227   if (useVertexBins)
1228     totalEvents = vertexAxis->GetNbins();
1229
1230   Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
1231   if (totalEvents > 0)
1232     totalTracks->Scale(1.0 / totalEvents);
1233   
1234   return totalTracks;
1235 }
1236
1237 //____________________________________________________________________
1238 TH1* AliUEHist::GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization)
1239 {
1240   // returns a histogram projected to pT,assoc with the provided cuts
1241   
1242    // unzoom all axes
1243   ResetBinLimits(fTrackHist[region]->GetGrid(step));
1244   ResetBinLimits(fEventHist->GetGrid(step));
1245   
1246   TH1D* tracks = 0;
1247   
1248     // the efficiency to have find an event depends on leading pT and this is not corrected for because anyway per bin we calculate tracks over events
1249   // therefore the number density must be calculated here per leading pT bin and then summed
1250
1251   if (multBinEnd > multBinBegin)
1252     Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
1253   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
1254   fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
1255   
1256   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
1257   fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
1258   
1259   // get real phi cuts due to binning
1260   Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
1261   Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
1262   Printf("phi min = %.2f and max = %.2f requested; due to binning range will be from %.2f to %.2f, i.e. a %.0f%% larger window", phiMin, phiMax, phiMinReal, phiMaxReal, (phiMaxReal - phiMinReal) / (phiMax - phiMin) * 100 - 100);
1263   
1264   Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
1265   Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
1266   
1267   TH1D* events = fEventHist->ShowProjection(0, step);
1268   
1269   for (Int_t bin=firstBin; bin<=lastBin; bin++)
1270   {
1271     fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
1272     
1273     // project to pT,assoc
1274     TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
1275     
1276     Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
1277     fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
1278     
1279     // normalize to get a density (deta dphi)
1280     Float_t normalization = 1;
1281     TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
1282     if (strcmp(axis->GetTitle(), "#eta") == 0)
1283     {
1284       Printf("Normalizing using eta-axis range");
1285       normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
1286     }
1287     else
1288     {
1289       Printf("Normalizating assuming accepted range of |eta| < 0.8");
1290       normalization *= 0.8 * 2;
1291     }
1292     
1293     // dphi
1294     if (!skipPhiNormalization)
1295       normalization *= phiMaxReal - phiMinReal;
1296     
1297     //Printf("Normalization: %f", normalization);
1298     tracksTmp->Scale(1.0 / normalization);
1299     
1300     // and now dpT (bins have different width)
1301     for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
1302     {
1303       tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1304       tracksTmp->SetBinError  (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1305     }
1306      
1307     Int_t nEvents = (Int_t) events->GetBinContent(bin);
1308     if (nEvents > 0)
1309       tracksTmp->Scale(1.0 / nEvents);
1310     Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
1311     
1312     if (!tracks)
1313       tracks = tracksTmp;
1314     else
1315     {
1316       tracks->Add(tracksTmp);
1317       delete tracksTmp;
1318     }
1319   }
1320   
1321   delete events;
1322
1323   ResetBinLimits(fTrackHist[region]->GetGrid(step));
1324   ResetBinLimits(fEventHist->GetGrid(step));
1325
1326   return tracks;
1327 }
1328
1329 void AliUEHist::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
1330 {
1331   // multiplies <grid> with <histogram> and put the result in <target>
1332   // <grid> has usually more dimensions than histogram. The axis which are used to choose the value 
1333   // from <histogram> are given in <var1> and <var2>
1334   //
1335   // if <histogram> is 0, just copies content from step1 to step2
1336   
1337   // clear target histogram
1338   if (grid != target)
1339     target->Reset();
1340   
1341   if (histogram != 0)
1342   {
1343     if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
1344       AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
1345       
1346     if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
1347       AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
1348   }
1349
1350   if (grid->GetNdimensions() > 6)
1351     AliFatal("Too many dimensions in THnSparse");
1352   
1353   Int_t bins[6];
1354   
1355   // optimized implementation
1356   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1357   {
1358     Double_t value = grid->GetBinContent(binIdx, bins);
1359     Double_t error = grid->GetBinError(binIdx);
1360     
1361     if (histogram != 0)
1362     {
1363       if (var2 < 0)
1364       {
1365         value *= histogram->GetBinContent(bins[var1]);
1366         error *= histogram->GetBinContent(bins[var1]);
1367       }
1368       else
1369       {
1370         value *= histogram->GetBinContent(bins[var1], bins[var2]);
1371         error *= histogram->GetBinContent(bins[var1], bins[var2]);
1372       }
1373     }
1374     
1375     target->SetBinContent(bins, value);
1376     target->SetBinError(bins, error);
1377   }
1378 }
1379
1380 //____________________________________________________________________
1381 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
1382 {
1383   // corrects from step1 to step2 by multiplying the tracks with trackCorrection
1384   // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
1385   // if var2 >= 0 a two dimension correction is assumed in trackCorrection
1386   //
1387   // if trackCorrection is 0, just copies content from step1 to step2
1388   
1389   for (UInt_t region=0; region<fkRegions; region++)
1390     CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
1391 }
1392
1393 //____________________________________________________________________
1394 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
1395 {
1396   //
1397   // see documentation of CorrectTracks above
1398   //
1399   
1400   if (!fTrackHist[region])
1401     return;
1402    
1403   THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
1404   THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
1405   
1406   Float_t entriesBefore = grid->GetEntries();
1407   
1408   MultiplyHistograms(grid, target, trackCorrection, var1, var2);
1409   
1410   Printf("AliUEHist::CorrectTracks: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0); 
1411 }
1412
1413 //____________________________________________________________________
1414 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
1415 {
1416   // corrects from step1 to step2 by multiplying the events with eventCorrection
1417   // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
1418   //
1419   // if eventCorrection is 0, just copies content from step1 to step2
1420   
1421   AliCFGridSparse* grid = fEventHist->GetGrid(step1);
1422   AliCFGridSparse* target = fEventHist->GetGrid(step2);
1423   
1424   Float_t entriesBefore = grid->GetEntries();
1425
1426   MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);
1427
1428   Printf("AliUEHist::CorrectEvents: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0); 
1429 }
1430
1431 //____________________________________________________________________
1432 void AliUEHist::Correct(AliUEHist* corrections)
1433 {
1434   // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
1435   //
1436   // in this object the data is expected in the step kCFStepReconstructed
1437   
1438   if (fHistogramType.Length() == 0)
1439   {
1440     Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
1441     if (fTrackHist[kToward]->GetNVar() < 5)
1442     {
1443       if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#varphid#eta") == 0)
1444         fHistogramType = "NumberDensitypT";
1445       else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#eta") == 0)
1446         fHistogramType = "SumpT";
1447     }
1448     else if (fTrackHist[kToward]->GetNVar() == 5)
1449     {
1450       if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
1451         fHistogramType = "NumberDensityPhi";
1452       else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
1453         fHistogramType = "NumberDensityPhiCentrality";
1454     }
1455       
1456     if (fHistogramType.Length() == 0)
1457       AliFatal("Cannot figure out histogram type. Try setting it manually...");
1458   }
1459   
1460   Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
1461   
1462   if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
1463   {
1464     // ---- track level
1465     
1466     // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
1467     // extracted as function of leading pT
1468     Bool_t biasFromMC = kFALSE;
1469     const char* projAxis = "z";
1470     Int_t secondBin = -1;
1471
1472     if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
1473     {
1474       projAxis = "yz";
1475       secondBin = 4;
1476     }
1477     
1478     for (UInt_t region = 0; region < fkRegions; region++)
1479     {
1480       if (!fTrackHist[region])
1481         continue;
1482    
1483       TH1* leadingBiasTracks =  0;
1484       if (biasFromMC)
1485       {
1486         leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
1487         Printf("WARNING: Using MC bias correction");
1488       }
1489       else
1490         leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1);          // from data
1491         
1492       CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
1493       if (region == kMin && fCombineMinMax)
1494       {
1495         CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
1496         delete leadingBiasTracks;
1497         break;
1498       }
1499       delete leadingBiasTracks;
1500     }
1501     
1502     TH1* leadingBiasEvents = 0;
1503     if (biasFromMC)
1504       leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
1505     else
1506       leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2);          // from data
1507       
1508     //new TCanvas; leadingBiasEvents->DrawCopy();
1509     CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
1510     delete leadingBiasEvents;
1511     
1512     // --- contamination correction ---
1513     
1514     TH2D* contamination = corrections->GetTrackingContamination();
1515     if (corrections->fContaminationEnhancement)
1516     {
1517       Printf("Applying contamination enhancement");
1518       
1519       for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
1520         for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
1521           contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
1522     }
1523     CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
1524     CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1525     delete contamination;
1526     
1527     // --- efficiency correction ---
1528     // correct single-particle efficiency for associated particles
1529     // in addition correct for efficiency on trigger particles (tracks AND events)
1530     
1531     TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
1532     // use kCFStepVertex as a temporary step (will be overwritten later anyway)
1533     CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1534     delete efficiencyCorrection;
1535     
1536     // correct pT,T
1537     efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
1538     CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
1539     CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
1540     delete efficiencyCorrection;
1541     
1542     // copy 
1543     CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
1544     CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
1545     
1546     // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
1547     // practically independent of low pT cut 
1548     TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
1549     
1550     // convert stage from true multiplicity to observed multiplicity by simple conversion factor
1551     TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
1552     vertexCorrectionObs->Reset();
1553     
1554     TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
1555     // some defaults
1556     func->SetParameters(0.1, 1, -0.7);
1557     vertexCorrection->Fit(func, "0I", "", 0, 3);
1558     for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
1559     {
1560       Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
1561       if (xPos < 3)
1562         vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
1563       else
1564         vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
1565     }
1566   
1567     #if 0
1568     new TCanvas;
1569     vertexCorrection->DrawCopy();
1570     vertexCorrectionObs->SetLineColor(2);
1571     vertexCorrectionObs->DrawCopy("same");
1572     func->SetRange(0, 4);
1573     func->DrawClone("same");
1574     #endif
1575     
1576     CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
1577     CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
1578     delete vertexCorrectionObs;
1579     delete vertexCorrection;
1580     delete func;
1581     
1582     // copy 
1583     CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
1584     CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
1585   }
1586   else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
1587   {
1588     if (fTrackHist[0]->GetNVar() <= 5)
1589     {
1590       // do corrections copying between steps
1591 //       CFStep step = kCFStepReconstructed;
1592       CFStep step = kCFStepBiasStudy;
1593       
1594       // copy 
1595       CorrectTracks(step, kCFStepTracked, 0, -1);
1596       CorrectEvents(step, kCFStepTracked, 0, -1);
1597       
1598       // Dont use eta in the following, because it is a Delta-eta axis
1599       
1600       // contamination correction
1601       // correct single-particle contamination for associated particles
1602       
1603       TH1* contamination = corrections->GetTrackingContamination(1);
1604       
1605       if (0)
1606       {
1607         Printf("Applying contamination enhancement");
1608         
1609         for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1610         {
1611           printf("%f", contamination->GetBinContent(bin));
1612           if (contamination->GetBinContent(bin) > 0)
1613             contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1614           printf(" --> %f\n", contamination->GetBinContent(bin));
1615         }
1616       }
1617         
1618       CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1619       delete contamination;    
1620       
1621       // correct for additional contamination due to trigger particle around phi ~ 0
1622       TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1623       if (0)
1624       {
1625         Printf("Applying contamination enhancement");
1626         
1627         for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1628           for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1629           {
1630             printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1631             if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1632               correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1633             printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1634           }
1635       }
1636       
1637 //       new TCanvas; correlatedContamination->DrawCopy("COLZ");
1638       CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1639 //       Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1640       
1641       delete correlatedContamination;
1642       
1643       // TODO correct for contamination of trigger particles (for tracks AND events)
1644       CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1645       
1646       // --- efficiency correction ---
1647       // correct single-particle efficiency for associated particles
1648       // in addition correct for efficiency on trigger particles (tracks AND events)
1649       
1650       // in bins of pT and centrality
1651       TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1652 //       new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1653       // use kCFStepAnaTopology as a temporary step 
1654       CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1655       delete efficiencyCorrection;
1656       
1657       // correct pT,T in bins of pT and centrality
1658       efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1659       CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1660       CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1661       delete efficiencyCorrection;
1662       
1663       // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1664       // copy 
1665       CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1666       CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1667     }
1668     else
1669     {
1670       // with 6 axes there is not enough memory, do the corrections in-place
1671       Printf(">>>>>>>> Applying corrections in place to reduce memory consumption");
1672       CFStep step = kCFStepBiasStudy;
1673 //       CFStep step = kCFStepReconstructed;
1674      
1675       // Dont use eta in the following, because it is a Delta-eta axis
1676       
1677       // --- contamination correction ---
1678       // correct single-particle contamination for associated particles
1679       // correct contamination for trigger particles (tracks AND events)
1680       
1681       // in bins of p,T
1682       TH1* contamination = corrections->GetTrackingContamination(1);
1683       
1684       if (0)
1685       {
1686         Printf("Applying contamination enhancement");
1687         
1688         for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1689         {
1690           printf("%f", contamination->GetBinContent(bin));
1691           if (contamination->GetBinContent(bin) > 0)
1692             contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1693           printf(" --> %f\n", contamination->GetBinContent(bin));
1694         }
1695       }
1696         
1697       // correct pT,A in bins of pT
1698       CorrectTracks(step, step, contamination, 1);
1699       delete contamination;
1700       
1701       // correct pT,T in bins of pT (for tracks AND events)
1702       contamination = corrections->GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, -1, 2);
1703       new TCanvas; contamination->DrawCopy();
1704       CorrectEvents(step, step, contamination, 0);
1705       CorrectTracks(step, step, contamination, 2);
1706       delete contamination;
1707       
1708       // correct for additional contamination due to trigger particle around phi ~ 0
1709       if (0)
1710       {
1711         TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1712         if (1)
1713         {
1714           Printf("Applying contamination enhancement");
1715           
1716           for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1717             for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1718             {
1719               printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1720               if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1721                 correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1722               printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1723             }
1724         }
1725
1726         // new TCanvas; correlatedContamination->DrawCopy("COLZ");
1727         CorrectCorrelatedContamination(step, 0, correlatedContamination);
1728
1729         delete correlatedContamination;
1730       }
1731       else
1732         Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1733       
1734       // --- tracking efficiency correction ---
1735       // correct single-particle efficiency for associated particles
1736       // correct for efficiency on trigger particles (tracks AND events)
1737       
1738       // in bins of pT and centrality
1739       TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1740       // new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1741
1742       // correct pT,A in bins of pT and centrality
1743       CorrectTracks(step, step, efficiencyCorrection, 1, 3);
1744       delete efficiencyCorrection;
1745       
1746       // correct pT,T in bins of pT and centrality
1747       efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1748       CorrectEvents(step, step, efficiencyCorrection, 0, 1);
1749       CorrectTracks(step, step, efficiencyCorrection, 2, 3);
1750       
1751       delete efficiencyCorrection;
1752     }
1753   }
1754   else
1755     AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1756 }
1757
1758 //____________________________________________________________________
1759 THnBase* AliUEHist::GetTrackEfficiencyND(CFStep step1, CFStep step2)
1760 {
1761   // creates a track-level efficiency by dividing step2 by step1
1762   // in all dimensions but the particle species one
1763   
1764   AliCFContainer* sourceContainer = fTrackHistEfficiency;
1765   // step offset because we start with kCFStepAnaTopology
1766   step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1767   step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1768   
1769   ResetBinLimits(sourceContainer->GetGrid(step1));
1770   ResetBinLimits(sourceContainer->GetGrid(step2));
1771
1772   if (fEtaMax > fEtaMin)
1773   {
1774     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1775     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1776     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1777   }
1778   
1779   Int_t dimensions[] = { 0, 1, 3, 4 };
1780   THnBase* generated = sourceContainer->GetGrid(step1)->GetGrid()->ProjectionND(4, dimensions);
1781   THnBase* measured = sourceContainer->GetGrid(step2)->GetGrid()->ProjectionND(4, dimensions);
1782   
1783 //   Printf("%d %d %f %f", step1, step2, generated->GetEntries(), measured->GetEntries());
1784   
1785   ResetBinLimits(sourceContainer->GetGrid(step1));
1786   ResetBinLimits(sourceContainer->GetGrid(step2));
1787
1788   measured->Divide(measured, generated, 1, 1, "B");
1789   
1790   delete generated;
1791   
1792   return measured;
1793 }
1794
1795 //____________________________________________________________________
1796 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1797 {
1798   // creates a track-level efficiency by dividing step2 by step1
1799   // projected to axis1 and axis2 (optional if >= 0)
1800   //
1801   // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1802   
1803   // integrate over regions
1804   // cache it for efficiency (usually more than one efficiency is requested)
1805   
1806   AliCFContainer* sourceContainer = 0;
1807   
1808   if (source == 0)
1809   {
1810     if (!fCache)
1811     {
1812       fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1813       for (UInt_t i = 1; i < fkRegions; i++)
1814         if (fTrackHist[i])
1815           fCache->Add(fTrackHist[i]);
1816     }
1817     sourceContainer = fCache;
1818   }
1819   else if (source == 1 || source == 2)
1820   {
1821     sourceContainer = fTrackHistEfficiency;
1822     // step offset because we start with kCFStepAnaTopology
1823     step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1824     step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1825   }
1826   else
1827     return 0;
1828         
1829   // reset all limits and set the right ones except those in axis1, axis2 and axis3
1830   ResetBinLimits(sourceContainer->GetGrid(step1));
1831   ResetBinLimits(sourceContainer->GetGrid(step2));
1832   if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1833   {
1834     Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1835     sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1836     sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1837   }
1838   if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1839   {
1840     Printf("Restricted pt-range to %f %f", fPtMin, fPtMax);
1841     sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1842     sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1843   }
1844   if (fPartSpecies != -1 && axis1 != 2 && axis2 != 2 && axis3 != 2)
1845   {
1846     Printf("Restricted to particle species %d", fPartSpecies);
1847     sourceContainer->GetGrid(step1)->SetRangeUser(2, fPartSpecies, fPartSpecies);
1848     sourceContainer->GetGrid(step2)->SetRangeUser(2, fPartSpecies, fPartSpecies);
1849   }
1850   if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1851   {
1852     Printf("Restricted centrality range to %f %f", fCentralityMin, fCentralityMax);
1853     sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1854     sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1855   }
1856   if (fZVtxMax > fZVtxMin && axis1 != 4 && axis2 != 4 && axis3 != 4)
1857   {
1858     Printf("Restricted z-vtx range to %f %f", fZVtxMin, fZVtxMax);
1859     sourceContainer->GetGrid(step1)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1860     sourceContainer->GetGrid(step2)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1861   }
1862   
1863   TH1* measured = 0;
1864   TH1* generated = 0;
1865     
1866   if (axis3 >= 0)
1867   {
1868     generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1869     measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1870   }
1871   else if (axis2 >= 0)
1872   {
1873     generated = sourceContainer->Project(step1, axis1, axis2);
1874     measured = sourceContainer->Project(step2, axis1, axis2);
1875   }
1876   else
1877   {
1878     generated = sourceContainer->Project(step1, axis1);
1879     measured = sourceContainer->Project(step2, axis1);
1880   }
1881   
1882   // check for bins with less than 50 entries, print warning
1883   Int_t binBegin[3];
1884   Int_t binEnd[3];
1885   
1886   binBegin[0] = 1;
1887   binBegin[1] = 1;
1888   binBegin[2] = 1;
1889   
1890   binEnd[0] = generated->GetNbinsX();
1891   binEnd[1] = generated->GetNbinsY();
1892   binEnd[2] = generated->GetNbinsZ();
1893   
1894   if (fEtaMax > fEtaMin)
1895   {
1896     if (axis1 == 0)
1897     {
1898       binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1899       binEnd[0]   = generated->GetXaxis()->FindBin(fEtaMax);
1900     }
1901     if (axis2 == 0)
1902     {
1903       binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1904       binEnd[1]   = generated->GetYaxis()->FindBin(fEtaMax);
1905     }
1906     if (axis3 == 0)
1907     {
1908       binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1909       binEnd[2]   = generated->GetZaxis()->FindBin(fEtaMax);
1910     }
1911   }
1912   
1913   if (fPtMax > fPtMin)
1914   {
1915     // TODO this is just checking up to 15 for now
1916     Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1917     if (axis1 == 1)
1918     {
1919       binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1920       binEnd[0]   = generated->GetXaxis()->FindBin(ptMax);
1921     }
1922     if (axis2 == 1)
1923     {
1924       binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1925       binEnd[1]   = generated->GetYaxis()->FindBin(ptMax);
1926     }
1927     if (axis3 == 1)
1928     {
1929       binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1930       binEnd[2]   = generated->GetZaxis()->FindBin(ptMax);
1931     }
1932   }
1933   
1934   Int_t total = 0;
1935   Int_t count = 0;
1936   Int_t vars[3];
1937   
1938   for (Int_t i=0; i<3; i++)
1939     vars[i] = binBegin[i];
1940     
1941   const Int_t limit = 50;
1942   while (1)
1943   {
1944     if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1945     {
1946       Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1947       count++;
1948     } 
1949     else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1950     {
1951       Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)", 
1952         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1953         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1954         generated->GetBinContent(vars[0], vars[1]));
1955       count++;
1956     }
1957     else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1958     {
1959       Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)", 
1960         generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1961         generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1962         generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1963         generated->GetBinContent(vars[0], vars[1], vars[2]));
1964       count++;
1965     }
1966     
1967     vars[2]++;
1968     if (vars[2] == binEnd[2]+1)
1969     {
1970       vars[2] = binBegin[2];
1971       vars[1]++;
1972     }
1973     
1974     if (vars[1] == binEnd[1]+1)
1975     {
1976       vars[1] = binBegin[1];
1977       vars[0]++;
1978     }
1979     
1980     if (vars[0] == binEnd[0]+1)
1981       break;
1982     total++;
1983   }
1984
1985   Printf("Correction has %d empty bins (out of %d bins)", count, total);
1986   
1987   // rebin if required
1988   if (source == 2)
1989   {
1990     TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1991     
1992     if (axis->GetNbins() < measured->GetNbinsX())
1993     {
1994       if (axis2 != -1)
1995       {
1996         // 2d rebin with variable axis does not exist in root
1997         
1998         TH1* tmp = measured;
1999         measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
2000         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
2001           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
2002           {
2003             ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
2004             measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
2005           }
2006         delete tmp;
2007         
2008         tmp = generated;
2009         generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
2010         for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
2011           for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
2012           {
2013             ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
2014             generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
2015           }
2016         delete tmp;
2017       }
2018       else
2019       {
2020         TH1* tmp = measured;
2021         measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2022         delete tmp;
2023         
2024         tmp = generated;
2025         generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2026         delete tmp;
2027       }
2028     }
2029     else if (axis->GetNbins() > measured->GetNbinsX())
2030     {
2031       if (axis2 != -1)
2032         AliFatal("Rebinning only works for 1d at present");
2033   
2034       // this is an unfortunate case where the number of bins has to be increased in principle
2035       // there is a region where the binning is finner in one histogram and a region where it is the other way round
2036       // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
2037       // only a certain binning is supported here
2038       if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
2039         AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
2040       
2041       Double_t newBins[] = {0.0, 0.5, 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};
2042       
2043       // reduce binning below 5 GeV/c
2044       TH1* tmp = measured;
2045       measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2046       delete tmp;
2047       
2048       // increase binning above 5 GeV/c
2049       tmp = measured;
2050       measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2051       for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
2052       {
2053         measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
2054         measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
2055       }
2056       delete tmp;
2057       
2058       // reduce binning below 5 GeV/c
2059       tmp = generated;
2060       generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2061       delete tmp;
2062       
2063       // increase binning above 5 GeV/c
2064       tmp = generated;
2065       generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2066       for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
2067       {
2068         generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
2069         generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
2070       }
2071       delete tmp;
2072     }
2073   }
2074   
2075   measured->Divide(measured, generated, 1, 1, "B");
2076   
2077   delete generated;
2078   
2079   ResetBinLimits(sourceContainer->GetGrid(step1));
2080   ResetBinLimits(sourceContainer->GetGrid(step2));
2081   
2082   return measured;
2083 }
2084
2085 //____________________________________________________________________
2086 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
2087 {
2088   // creates a event-level efficiency by dividing step2 by step1
2089   // projected to axis1 and axis2 (optional if >= 0)
2090   
2091   if (ptLeadMax > ptLeadMin)
2092   {
2093     fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2094     fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2095   }
2096   
2097   TH1* measured = 0;
2098   TH1* generated = 0;
2099     
2100   if (axis2 >= 0)
2101   {
2102     generated = fEventHist->Project(step1, axis1, axis2);
2103     measured = fEventHist->Project(step2, axis1, axis2);
2104   }
2105   else
2106   {
2107     generated = fEventHist->Project(step1, axis1);
2108     measured = fEventHist->Project(step2, axis1);
2109   }
2110   
2111   measured->Divide(measured, generated, 1, 1, "B");
2112   
2113   delete generated;
2114   
2115   if (ptLeadMax > ptLeadMin)
2116   {
2117     fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
2118     fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
2119   }
2120   
2121   return measured;
2122 }
2123
2124 //____________________________________________________________________
2125 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
2126 {
2127   // weights each entry of the 3d histogram hist1 with the 1d histogram hist2 
2128   // where the matching is done of the z axis of hist1 with the x axis of hist2
2129   
2130   if (hist1->GetNbinsZ() != hist2->GetNbinsX())
2131     AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
2132   
2133   for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
2134   {
2135     for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
2136     {
2137       for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
2138       {
2139         if (hist2->GetBinContent(z) > 0)
2140         {
2141           hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
2142           hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
2143         }
2144         else
2145         {
2146           hist1->SetBinContent(x, y, z, 0);
2147           hist1->SetBinError(x, y, z, 0);
2148         }
2149       }
2150     }
2151   }
2152 }  
2153
2154 //____________________________________________________________________
2155 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
2156 {
2157   // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
2158   // in the given region (sum over all regions is calculated if region == -1)
2159   // done by weighting the track-level distribution with the number of events as function of leading pT
2160   // and then calculating the ratio between the distributions
2161   // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
2162   //   no projection is done if axis == 0
2163   // weighting: 0 = tracks weighted with events (as discussed above)
2164   //            1 = only track bias is returned
2165   //            2 = only event bias is returned
2166   
2167   AliCFContainer* tmp = 0;
2168   
2169   if (region == -1)
2170   {
2171     tmp = (AliCFContainer*) fTrackHist[0]->Clone();
2172     for (UInt_t i = 1; i < fkRegions; i++)
2173       if (fTrackHist[i])
2174         tmp->Add(fTrackHist[i]);
2175   }
2176   else if (region == kMin && fCombineMinMax)
2177   {
2178     tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
2179     tmp->Add(fTrackHist[kMax]);
2180   }
2181   else
2182     tmp = fTrackHist[region];
2183   
2184   ResetBinLimits(tmp->GetGrid(step1));
2185   ResetBinLimits(fEventHist->GetGrid(step1));
2186   SetBinLimits(tmp->GetGrid(step1));
2187   
2188   ResetBinLimits(tmp->GetGrid(step2));
2189   ResetBinLimits(fEventHist->GetGrid(step2));
2190   SetBinLimits(tmp->GetGrid(step2));
2191   
2192   TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
2193   TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
2194   if (weighting == 0)
2195     WeightHistogram(hist1, events1);
2196   
2197   TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
2198   TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
2199   if (weighting == 0)
2200     WeightHistogram(hist2, events2);
2201   
2202   TH1* generated = hist1;
2203   TH1* measured = hist2;
2204   
2205   if (weighting == 0 || weighting == 1)
2206   {
2207     if (axis)
2208     {
2209       if (leadPtMax > leadPtMin)
2210       {
2211         hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2212         hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2213       }
2214       
2215       if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
2216       {
2217         hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2218         hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2219       }
2220     
2221       generated = hist1->Project3D(axis);
2222       measured  = hist2->Project3D(axis);
2223       
2224       // delete hists here if projection has been done
2225       delete hist1;
2226       delete hist2;
2227     }
2228     delete events1;
2229     delete events2;
2230   }
2231   else if (weighting == 2)
2232   {
2233     delete hist1;
2234     delete hist2;
2235     generated = events1;
2236     measured = events2;
2237   }
2238   
2239   measured->Divide(generated);
2240   
2241   delete generated;
2242   
2243   ResetBinLimits(tmp->GetGrid(step1));
2244   ResetBinLimits(tmp->GetGrid(step2));
2245   
2246   if ((region == -1) || (region == kMin && fCombineMinMax))
2247     delete tmp;
2248   
2249   return measured;
2250 }
2251
2252 //____________________________________________________________________
2253 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
2254 {
2255   // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
2256   
2257   if (!fTrackHist[region])
2258     return;
2259    
2260   THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
2261   
2262   Int_t var1 = 1;
2263   Int_t var2 = 2;
2264   
2265   if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
2266     AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
2267     
2268   if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
2269     AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
2270   
2271   // optimized implementation
2272   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2273   {
2274     Int_t bins[6];
2275     
2276     Double_t value = grid->GetBinContent(binIdx, bins);
2277     Double_t error = grid->GetBinError(binIdx);
2278     
2279     // check eta and phi axes
2280     if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
2281       continue;
2282     if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
2283       continue;
2284     
2285     value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2286     error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2287     
2288     grid->SetBinContent(bins, value);
2289     grid->SetBinError(bins, error);
2290   }
2291  
2292   Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
2293 }
2294
2295 //____________________________________________________________________
2296 TH2* AliUEHist::GetCorrelatedContamination()
2297 {
2298   // contamination correlated with the trigger particle is evaluated between step kCFStepTracked and kCFStepTrackedOnlyPrim in the region of delta eta and delta phi < 0.1 (smallest bin!)
2299   
2300   Int_t step1 = kCFStepTrackedOnlyPrim;
2301   Int_t step2 = kCFStepTracked;
2302   
2303   fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
2304   fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
2305   TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
2306   
2307   fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
2308   fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
2309   TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
2310   
2311   tracksStep1->Divide(tracksStep2);
2312   
2313   TH1* triggersStep1 = fEventHist->Project(step1, 0);
2314   TH1* triggersStep2 = fEventHist->Project(step2, 0);
2315   
2316   TH1* singleParticle = GetTrackingContamination(1);
2317   
2318   for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
2319     for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
2320       if (singleParticle->GetBinContent(x) > 0 && triggersStep1->GetBinContent(y) > 0)
2321         tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
2322       else
2323         tracksStep1->SetBinContent(x, y, 0);
2324         
2325   delete singleParticle;
2326   delete tracksStep2;
2327   delete triggersStep1;
2328   delete triggersStep2;
2329         
2330   return tracksStep1;
2331 }
2332
2333 //____________________________________________________________________
2334 TH2D* AliUEHist::GetTrackingEfficiency()
2335 {
2336   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2337   // integrates over the regions and all other variables than pT and eta to increase the statistics
2338   //
2339   // returned histogram has to be deleted by the user
2340
2341   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
2342 }
2343   
2344 //____________________________________________________________________
2345 TH2D* AliUEHist::GetFakeRate()
2346 {
2347   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), 0, 1));
2348 }
2349
2350 //____________________________________________________________________
2351 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
2352 {
2353   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2354   // integrates over the regions and all other variables than pT, centrality to increase the statistics
2355   //
2356   // returned histogram has to be deleted by the user
2357
2358   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
2359 }
2360
2361 //____________________________________________________________________
2362 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
2363 {
2364   // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2365   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2366
2367   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
2368 }
2369
2370 //____________________________________________________________________
2371 TH1D* AliUEHist::GetFakeRate(Int_t axis)
2372 {
2373   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), axis));
2374 }
2375 //____________________________________________________________________
2376 TH2D* AliUEHist::GetTrackingCorrection()
2377 {
2378   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2379   // integrates over the regions and all other variables than pT and eta to increase the statistics
2380   //
2381   // returned histogram has to be deleted by the user
2382
2383   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
2384 }
2385   
2386 //____________________________________________________________________
2387 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
2388 {
2389   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2390   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2391
2392   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
2393 }
2394
2395 //____________________________________________________________________
2396 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
2397 {
2398   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2399   // integrates over the regions and all other variables than pT and eta to increase the statistics
2400   //
2401   // returned histogram has to be deleted by the user
2402
2403   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
2404 }
2405   
2406 //____________________________________________________________________
2407 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
2408 {
2409   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2410   // integrates over the regions and all other variables than pT and centrality to increase the statistics
2411   //
2412   // returned histogram has to be deleted by the user
2413
2414   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
2415 }
2416   
2417 //____________________________________________________________________
2418 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
2419 {
2420   // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2421   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2422
2423   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
2424 }
2425
2426 //____________________________________________________________________
2427 TH2D* AliUEHist::GetTrackingContamination()
2428 {
2429   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2430   // integrates over the regions and all other variables than pT and eta to increase the statistics
2431   //
2432   // returned histogram has to be deleted by the user
2433
2434   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
2435 }
2436   
2437 //____________________________________________________________________
2438 TH2D* AliUEHist::GetTrackingContaminationCentrality()
2439 {
2440   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2441   // integrates over the regions and all other variables than pT and centrality to increase the statistics
2442   //
2443   // returned histogram has to be deleted by the user
2444
2445   return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
2446 }
2447   
2448 //____________________________________________________________________
2449 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
2450 {
2451   // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2452   // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2453
2454   return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
2455 }
2456
2457 //____________________________________________________________________
2458 const char* AliUEHist::GetRegionTitle(Region region)
2459 {
2460   // returns the name of the given region
2461   
2462   switch (region)
2463   {
2464     case kToward:
2465       return "Towards";
2466     case kAway:
2467       return "Away";
2468     case kMin:
2469       return (fCombineMinMax) ? "Transverse" : "Min";
2470     case kMax:
2471       return "Max";
2472   }
2473   
2474   return 0;
2475 }
2476   
2477 //____________________________________________________________________
2478 const char* AliUEHist::GetStepTitle(CFStep step)
2479 {
2480   // returns the name of the given step
2481   
2482   switch (step)
2483   {
2484     case kCFStepAll:
2485       return "All events";
2486     case kCFStepTriggered:
2487       return "Triggered";
2488     case kCFStepVertex:
2489       return "Primary Vertex";
2490     case kCFStepAnaTopology:
2491       return "Required analysis topology";
2492     case kCFStepTrackedOnlyPrim:
2493       return "Tracked (matched MC, only primaries)";
2494     case kCFStepTracked:
2495       return "Tracked (matched MC, all)";
2496     case kCFStepReconstructed:
2497       return "Reconstructed";
2498     case kCFStepRealLeading:
2499       return "Correct leading particle identified";
2500     case kCFStepBiasStudy:
2501       return "Bias study applying tracking efficiency";
2502     case kCFStepBiasStudy2:
2503       return "Bias study applying tracking efficiency in two steps";
2504     case kCFStepCorrected:
2505       return "Corrected for efficiency on-the-fly";
2506   }
2507   
2508   return 0;
2509 }
2510
2511 //____________________________________________________________________
2512 void AliUEHist::CopyReconstructedData(AliUEHist* from)
2513 {
2514   // copies those histograms extracted from ESD to this object
2515   
2516   // TODO at present only the pointers are copied
2517   
2518   for (Int_t region=0; region<4; region++)
2519   {
2520     if (!fTrackHist[region])
2521       continue;
2522   
2523     fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
2524     //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2525     fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
2526   }
2527     
2528   fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
2529   //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2530   fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy,     from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
2531 }
2532
2533 void AliUEHist::DeepCopy(AliUEHist* from)
2534 {
2535   // copies the entries of this object's members from the object <from> to this object
2536   // fills using the fill function and thus allows that the objects have different binning
2537
2538   for (Int_t region=0; region<4; region++)
2539   {
2540     if (!fTrackHist[region] || !from->fTrackHist[region])
2541       continue;
2542   
2543     for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
2544     {
2545       Printf("Copying region %d step %d", region, step);
2546       THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
2547       THnSparse* source = from->fTrackHist[region]->GetGrid(step)->GetGrid();
2548       
2549       target->Reset();
2550       target->RebinnedAdd(source);
2551     }
2552   }
2553   
2554   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2555   {
2556     Printf("Ev: Copying step %d", step);
2557     THnSparse* target = fEventHist->GetGrid(step)->GetGrid();
2558     THnSparse* source = from->fEventHist->GetGrid(step)->GetGrid();
2559
2560     target->Reset();
2561     target->RebinnedAdd(source);
2562   }
2563   
2564   for (Int_t step=0; step<TMath::Min(fTrackHistEfficiency->GetNStep(), from->fTrackHistEfficiency->GetNStep()); step++)
2565   {
2566     if (!from->fTrackHistEfficiency->GetGrid(step))
2567       continue;
2568     
2569     Printf("Eff: Copying step %d", step);
2570     THnSparse* target = fTrackHistEfficiency->GetGrid(step)->GetGrid();
2571     THnSparse* source = from->fTrackHistEfficiency->GetGrid(step)->GetGrid();
2572
2573     target->Reset();
2574     target->RebinnedAdd(source);
2575   }
2576 }
2577
2578 void AliUEHist::SymmetrizepTBins()
2579 {
2580   // copy pt,a < pt,t bins to pt,a > pt,t (inverting deltaphi and delta eta as it should be) including symmetric bins
2581   
2582   for (Int_t region=0; region<4; region++)
2583   {
2584     if (!fTrackHist[region])
2585       continue;
2586   
2587     for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
2588     {
2589       Printf("Copying region %d step %d", region, step);
2590       THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
2591       if (target->GetEntries() == 0)
2592         continue;
2593       
2594       // for symmetric bins
2595       THnSparse* source = (THnSparse*) target->Clone();
2596       
2597       Int_t zVtxBins = 1;
2598       if (target->GetNdimensions() > 5)
2599         zVtxBins = target->GetAxis(5)->GetNbins();
2600       
2601       // axes: 0 delta eta; 1 pT,a; 2 pT,t; 3 centrality; 4 delta phi; 5 vtx-z
2602       for (Int_t i3 = 1; i3 <= target->GetAxis(3)->GetNbins(); i3++)
2603         for (Int_t i5 = 1; i5 <= zVtxBins; i5++)
2604         {
2605           for (Int_t i1 = 1; i1 <= target->GetAxis(1)->GetNbins(); i1++)
2606             for (Int_t i2 = 1; i2 <= target->GetAxis(2)->GetNbins(); i2++)
2607             {
2608               // find source bin
2609               Int_t binA = target->GetAxis(1)->FindBin(target->GetAxis(2)->GetBinCenter(i2));
2610               Int_t binT = target->GetAxis(2)->FindBin(target->GetAxis(1)->GetBinCenter(i1));
2611               
2612               Printf("(%d %d) Copying from %d %d to %d %d", i3, i5, binA, binT, i1, i2);
2613               
2614               for (Int_t i0 = 1; i0 <= target->GetAxis(0)->GetNbins(); i0++)
2615                 for (Int_t i4 = 1; i4 <= target->GetAxis(4)->GetNbins(); i4++)
2616                 {
2617                   Int_t binEta = target->GetAxis(0)->FindBin(-target->GetAxis(0)->GetBinCenter(i0));
2618                   Double_t phi = -target->GetAxis(4)->GetBinCenter(i4);
2619                   if (phi < -TMath::Pi()/2)
2620                     phi += TMath::TwoPi();
2621                   Int_t binPhi = target->GetAxis(4)->FindBin(phi);
2622                   
2623                   Int_t binSource[] = { binEta, binA, binT, i3, binPhi, i5 };
2624                   Int_t binTarget[] = { i0, i1, i2, i3, i4, i5 };
2625                   
2626                   Double_t value = source->GetBinContent(binSource);
2627                   Double_t error = source->GetBinError(binSource);
2628                   
2629                   if (error == 0)
2630                     continue;
2631                   
2632                   Double_t value2 = target->GetBinContent(binTarget);
2633                   Double_t error2 = target->GetBinError(binTarget);
2634                   
2635                   Double_t sum = value;
2636                   Double_t err = error;
2637                   
2638                   if (error2 > 0)
2639                   {
2640                     sum = value + value2;
2641                     err = TMath::Sqrt(error * error + error2 * error2);
2642                   }
2643
2644                   // Printf("  Values: %f +- %f; %f +- %f --> %f +- %f", value, error, value2, error2, sum, err);
2645                   
2646                   target->SetBinContent(binTarget, sum);
2647                   target->SetBinError(binTarget, err);
2648                 }
2649             }
2650         }
2651         
2652         delete source;
2653     }
2654   }
2655 }
2656
2657 //____________________________________________________________________
2658 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
2659 {
2660   // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
2661
2662   Float_t fitRangeBegin = 5.01;
2663   Float_t fitRangeEnd = 14.99;
2664   Float_t extendRangeBegin = 10.01;
2665
2666   if (fTrackHistEfficiency->GetNVar() == 3)
2667   {
2668     TH1* obj = GetTrackingEfficiency(1);
2669   
2670     if (verbose)
2671     {
2672       new TCanvas; 
2673       obj->Draw();
2674     }
2675     
2676     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2677   
2678     Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
2679   
2680     obj = GetTrackingContamination(1);
2681   
2682     if (verbose)
2683     {
2684       new TCanvas; 
2685       obj->Draw();
2686     }
2687     
2688     obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2689   
2690     Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
2691   
2692     Printf("AliUEHist::ExtendTrackingEfficiency: Fitted efficiency between %f and %f and got %f tracking efficiency and %f tracking contamination correction. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, trackingEff, trackingCont, extendRangeBegin, fEtaMin, fEtaMax);
2693   
2694     // extend for full pT range
2695     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2696       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2697         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2698         {
2699           
2700           Int_t bins[3];
2701           bins[0] = x;
2702           bins[1] = y;
2703           bins[2] = z;
2704           
2705           fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2706           fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
2707           fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
2708         }
2709   }
2710   else if (fTrackHistEfficiency->GetNVar() == 4)
2711   {
2712     // fit in centrality intervals of 20% for efficiency, one bin for contamination
2713     Float_t* trackingEff = 0;
2714     Float_t* trackingCont = 0;
2715     Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
2716     Int_t nCentralityBins = 5;
2717     
2718     Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
2719     
2720     // 0 = eff; 1 = cont
2721     for (Int_t caseNo = 0; caseNo < 2; caseNo++)
2722     {
2723       Float_t* target = 0;
2724       Int_t centralityBinsLocal = nCentralityBins;
2725       
2726       if (caseNo == 0)
2727       {
2728         trackingEff = new Float_t[centralityBinsLocal];
2729         target = trackingEff;
2730       }
2731       else
2732       {
2733         centralityBinsLocal = 1;
2734         trackingCont = new Float_t[centralityBinsLocal];
2735         target = trackingCont;
2736       }
2737     
2738       for (Int_t i=0; i<centralityBinsLocal; i++)
2739       {
2740         if (centralityBinsLocal == 1)
2741           SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
2742         else
2743           SetCentralityRange(centralityBins[i] + 0.1, centralityBins[i+1] - 0.1);
2744         TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
2745         if (verbose)
2746         {
2747           new TCanvas;
2748           proj->DrawCopy();
2749         }
2750         if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
2751           target[i] = proj->GetFunction("pol0")->GetParameter(0);
2752         else
2753           target[i] = 0;
2754         Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
2755       }
2756     }
2757   
2758     // extend for full pT range
2759     for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2760       for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2761         for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2762         {
2763           for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
2764           {
2765             
2766             Int_t bins[4];
2767             bins[0] = x;
2768             bins[1] = y;
2769             bins[2] = z;
2770             bins[3] = z2;
2771             
2772             Int_t z2Bin = 0;
2773             while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
2774               z2Bin++;
2775             
2776             //Printf("%d %d", z2, z2Bin);
2777             
2778             fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2779             fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
2780             if (trackingCont[0] > 0)
2781               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
2782             else  
2783               fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
2784           }
2785        }
2786        
2787      delete[] trackingEff;
2788      delete[] trackingCont;
2789    }
2790    
2791    SetCentralityRange(0, 0);
2792 }
2793
2794 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
2795 {
2796   // corrects the dphi distribution with an extra factor close to dphi ~ 0
2797
2798   Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
2799
2800   THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
2801   
2802   // optimized implementation
2803   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2804   {
2805     Int_t bins[6];
2806     Double_t value = grid->GetBinContent(binIdx, bins);
2807     Double_t error = grid->GetBinError(binIdx);
2808     
2809     Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
2810     if (TMath::Abs(binCenter) < 0.2)
2811     {
2812       value *= 0.985;
2813       error *= 0.985;
2814     }
2815     else if (TMath::Abs(binCenter) < 0.3)
2816     {
2817       value *= 0.9925;
2818       error *= 0.9925;
2819     }
2820     
2821     grid->SetBinContent(bins, value);
2822     grid->SetBinError(bins, error);
2823   }
2824 }
2825
2826 void AliUEHist::Scale(Double_t factor)
2827 {
2828   // scales all contained histograms by the given factor
2829   
2830   for (Int_t i=0; i<4; i++)
2831     if (fTrackHist[i])
2832       fTrackHist[i]->Scale(factor);
2833   
2834   fEventHist->Scale(factor);
2835   fTrackHistEfficiency->Scale(factor);
2836 }
2837
2838 void AliUEHist::Reset()
2839 {
2840   // resets all contained histograms
2841   
2842   for (Int_t i=0; i<4; i++)
2843     if (fTrackHist[i])
2844       for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
2845         fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2846   
2847   for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2848     fEventHist->GetGrid(step)->GetGrid()->Reset();
2849     
2850   for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
2851     fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
2852 }
2853
2854 THnBase* AliUEHist::ChangeToThn(THnBase* sparse)
2855 {
2856   // change the object to THn for faster processing
2857   
2858         // convert to THn (SEGV's for some strange reason...) 
2859         // x = THn::CreateHn("a", "a", sparse);
2860         
2861   // own implementation
2862   Int_t nBins[10];
2863   for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2864     nBins[i] = sparse->GetAxis(i)->GetNbins();
2865   THn* tmpTHn = new THnF(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse->GetNdimensions(), nBins, 0, 0);
2866   for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2867   {
2868     tmpTHn->SetBinEdges(i, sparse->GetAxis(i)->GetXbins()->GetArray());
2869     tmpTHn->GetAxis(i)->SetTitle(sparse->GetAxis(i)->GetTitle());
2870   }
2871   tmpTHn->RebinnedAdd(sparse);
2872   
2873   return tmpTHn;
2874 }
2875
2876 void AliUEHist::CondenseBin(THnSparse* grid, THnSparse* target, Int_t axis, Float_t targetValue, Float_t from, Float_t to)
2877 {
2878   //
2879   // loops through the histogram and moves all entries to a single point <targetValue> on the axis <axis>
2880   // if <from> and <to> are set, then moving only occurs if value on <axis> is betweem <from> and <to>
2881   //
2882
2883   if (grid->GetNdimensions() > 6)
2884     AliFatal("Too many dimensions in THnSparse");
2885   
2886   Int_t targetBin = grid->GetAxis(axis)->FindBin(targetValue);
2887   AliInfo(Form("Target bin on axis %d with value %f is %d", axis, targetValue, targetBin));
2888   
2889   Int_t fromBin = 1;
2890   Int_t toBin = grid->GetAxis(axis)->GetNbins();
2891   if (to > from)
2892   {
2893     fromBin = grid->GetAxis(axis)->FindBin(from);
2894     toBin = grid->GetAxis(axis)->FindBin(to);
2895     AliInfo(Form("Only condensing between bin %d and %d", fromBin, toBin));
2896   }
2897   
2898   Int_t bins[6];
2899   for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2900   {
2901     Double_t value = grid->GetBinContent(binIdx, bins);
2902     Double_t error = grid->GetBinError(binIdx);
2903     
2904     if (bins[axis] >= fromBin && bins[axis] <= toBin)
2905       bins[axis] = targetBin;
2906
2907     value += target->GetBinContent(bins);
2908     error = TMath::Sqrt(error * error + target->GetBinError(bins) * target->GetBinError(bins));
2909     
2910     target->SetBinContent(bins, value);
2911     target->SetBinError(bins, error);
2912   }
2913 }
2914
2915 void AliUEHist::CondenseBin(CFStep step, Int_t trackAxis, Int_t eventAxis, Float_t targetValue, Float_t from, Float_t to, CFStep tmpStep)
2916 {
2917   // loops through the histogram at <step> and moves all entries to a single point <targetValue> on the axes 
2918   // <trackAxis> and <eventAxis>. <tmpStep> is used to temporary store the data
2919   // This is useful e.g. to move bin content around for MC productions where the centrality selection did
2920   // not yield the desired result
2921
2922   // reset tmpStep
2923   fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2924   for (UInt_t i=0; i<fkRegions; i++)
2925     if (fTrackHist[i])
2926       fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2927
2928   // copy to tmpStep
2929   CorrectTracks(step, tmpStep, 0, -1);
2930   CorrectEvents(step, tmpStep, 0, -1);
2931
2932   // reset step
2933   fEventHist->GetGrid(step)->GetGrid()->Reset();
2934   for (UInt_t i=0; i<fkRegions; i++)
2935     if (fTrackHist[i])
2936       fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2937   
2938   // rewriting
2939   for (UInt_t i=0; i<fkRegions; i++)
2940   {
2941     if (!fTrackHist[i])
2942       continue;
2943     
2944     THnSparse* grid = fTrackHist[i]->GetGrid(tmpStep)->GetGrid();
2945     THnSparse* target = fTrackHist[i]->GetGrid(step)->GetGrid();
2946     
2947     CondenseBin(grid, target, trackAxis, targetValue, from, to);
2948   }
2949   CondenseBin(fEventHist->GetGrid(tmpStep)->GetGrid(), fEventHist->GetGrid(step)->GetGrid(), eventAxis, targetValue, from, to);
2950   
2951   // reset tmpStep
2952   fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2953   for (UInt_t i=0; i<fkRegions; i++)
2954     if (fTrackHist[i])
2955       fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2956 }