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