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