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