]>
Commit | Line | Data |
---|---|---|
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 | |
40 | ClassImp(AliUEHist) | |
41 | ||
408d1ac9 | 42 | const Int_t AliUEHist::fgkCFSteps = 11; |
a75aacd6 | 43 | |
3f3f12d9 | 44 | AliUEHist::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 | |
340 | TString 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 | 363 | Double_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 | //_____________________________________________________________________________ |
400 | AliUEHist::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 | //____________________________________________________________________ |
434 | void 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 | //____________________________________________________________________ | |
443 | AliUEHist::~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 | //____________________________________________________________________ | |
482 | AliUEHist &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 | //____________________________________________________________________ | |
493 | void 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 | //____________________________________________________________________ | |
534 | Long64_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 | 592 | void 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 | //____________________________________________________________________ |
605 | void 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 | //____________________________________________________________________ | |
613 | void 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 | //____________________________________________________________________ |
623 | void AliUEHist::ResetBinLimits(AliCFGridSparse* grid) | |
624 | { | |
625 | // resets all bin limits | |
626 | ||
d4b3dbfc | 627 | ResetBinLimits(grid->GetGrid()); |
a75aacd6 | 628 | } |
d4b3dbfc | 629 | |
144bd037 | 630 | //____________________________________________________________________ |
631 | void 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 | 708 | TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm, Long64_t* normEvents) |
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 | //____________________________________________________________________ |
860 | void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist) | |
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 | //____________________________________________________________________ |
895 | void 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 | 949 | TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t normalizePerTrigger, 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 | 1251 | TH1* 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 | /* |
1276 | TH2* AliUEHist::GetPtDistInPhiRegion(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiBegin, Float_t phiEnd) | |
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 | //____________________________________________________________________ |
1287 | TH2* AliUEHist::GetSumOfRatios(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm, Bool_t useVertexBins) | |
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 | //____________________________________________________________________ |
1390 | TH1* AliUEHist::GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization) | |
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 | 1481 | void 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 | //____________________________________________________________________ | |
1533 | void 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 | //____________________________________________________________________ | |
1546 | void 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 | 1566 | void 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 | //____________________________________________________________________ | |
1584 | void 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 | //____________________________________________________________________ |
1911 | THnBase* 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 | 1948 | TH1* 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 | //____________________________________________________________________ | |
2238 | TH1* 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 | //____________________________________________________________________ | |
2277 | void 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 | 2307 | TH1* 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 | //____________________________________________________________________ |
2405 | void 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 | //____________________________________________________________________ | |
2448 | TH2* 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 | //____________________________________________________________________ |
2486 | TH2D* 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 | //____________________________________________________________________ |
2497 | TH2D* AliUEHist::GetFakeRate() | |
2498 | { | |
59375a30 | 2499 | return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), 0, 1)); |
b591fb9c | 2500 | } |
2501 | ||
2a910c25 | 2502 | //____________________________________________________________________ |
2503 | TH2D* 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 | //____________________________________________________________________ |
2514 | TH1D* 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 | //____________________________________________________________________ |
2523 | TH1D* AliUEHist::GetFakeRate(Int_t axis) | |
2524 | { | |
59375a30 | 2525 | return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, (CFStep) (kCFStepTracked+3), axis)); |
46848f0b | 2526 | } |
a75aacd6 | 2527 | //____________________________________________________________________ |
2528 | TH2D* 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 | //____________________________________________________________________ | |
2539 | TH1D* 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 | //____________________________________________________________________ |
2548 | TH2D* 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 | //____________________________________________________________________ |
2559 | TH2D* 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 | //____________________________________________________________________ |
2570 | TH1D* 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 | //____________________________________________________________________ |
2579 | TH2D* 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 | //____________________________________________________________________ |
2590 | TH2D* 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 | //____________________________________________________________________ |
2601 | TH1D* 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 | //____________________________________________________________________ | |
2610 | const 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 | //____________________________________________________________________ | |
2630 | const 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 | //____________________________________________________________________ | |
2664 | void 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 | 2685 | void 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 | 2730 | void 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 | 2810 | void 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 | |
2946 | void 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 | |
2978 | void 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 | ||
2990 | void 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 | |
3006 | THnBase* 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 | 3028 | void 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 | 3067 | void 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 | } |