]>
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" |
a75aacd6 | 37 | |
38 | ClassImp(AliUEHist) | |
39 | ||
40 | const Int_t AliUEHist::fgkCFSteps = 10; | |
41 | ||
42 | AliUEHist::AliUEHist(const char* reqHist) : | |
d1c75d06 | 43 | TObject(), |
a75aacd6 | 44 | fkRegions(4), |
45 | fEventHist(0), | |
b1831bcb | 46 | fTrackHistEfficiency(0), |
a75aacd6 | 47 | fEtaMin(0), |
48 | fEtaMax(0), | |
49 | fPtMin(0), | |
50 | fPtMax(0), | |
144bd037 | 51 | fContaminationEnhancement(0), |
a75aacd6 | 52 | fCombineMinMax(0), |
53 | fCache(0) | |
54 | { | |
55 | // Constructor | |
56 | ||
57 | for (Int_t i=0; i<fkRegions; i++) | |
58 | fTrackHist[i] = 0; | |
59 | ||
60 | if (strlen(reqHist) == 0) | |
61 | return; | |
62 | ||
df269636 | 63 | AliLog::SetClassDebugLevel("AliCFContainer", -1); |
e741fadd | 64 | AliLog::SetClassDebugLevel("AliCFGridSparse", -3); |
df269636 | 65 | |
a75aacd6 | 66 | const char* title = ""; |
67 | ||
68 | // track level | |
69 | Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity | |
70 | Int_t iTrackBin[5]; | |
71 | Double_t* trackBins[5]; | |
72 | const char* trackAxisTitle[5]; | |
73 | ||
74 | // eta | |
75 | iTrackBin[0] = 20; | |
76 | Double_t etaBins[20+1]; | |
77 | for (Int_t i=0; i<=iTrackBin[0]; i++) | |
78 | etaBins[i] = -1.0 + 0.1 * i; | |
79 | trackBins[0] = etaBins; | |
80 | trackAxisTitle[0] = "#eta"; | |
81 | ||
e0331fd9 | 82 | // delta eta |
83 | const Int_t kNDeltaEtaBins = 20; | |
84 | Double_t deltaEtaBins[kNDeltaEtaBins+1]; | |
85 | for (Int_t i=0; i<=kNDeltaEtaBins; i++) | |
86 | deltaEtaBins[i] = -2.0 + 0.2 * i; | |
87 | ||
a75aacd6 | 88 | // pT |
89 | iTrackBin[1] = 39; | |
90 | Double_t pTBins[] = {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 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}; | |
91 | trackBins[1] = pTBins; | |
92 | trackAxisTitle[1] = "p_{T} (GeV/c)"; | |
93 | ||
94 | // pT,lead binning 1 | |
95 | const Int_t kNLeadingpTBins = 100; | |
96 | Double_t leadingpTBins[kNLeadingpTBins+1]; | |
97 | for (Int_t i=0; i<=kNLeadingpTBins; i++) | |
98 | leadingpTBins[i] = 0.5 * i; | |
99 | ||
100 | // pT,lead binning 2 | |
ada1a03f | 101 | const Int_t kNLeadingpTBins2 = 13; |
102 | Double_t leadingpTBins2[] = { 0.0, 0.5, 1.0, 2.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 100.0 }; | |
a75aacd6 | 103 | |
104 | // phi,lead | |
105 | const Int_t kNLeadingPhiBins = 40; | |
106 | Double_t leadingPhiBins[kNLeadingPhiBins+1]; | |
107 | for (Int_t i=0; i<=kNLeadingPhiBins; i++) | |
2ac8dc5c | 108 | leadingPhiBins[i] = -0.5 * TMath::Pi() + 1.0 / 40 * i * TMath::TwoPi(); |
a75aacd6 | 109 | |
110 | // multiplicity | |
111 | const Int_t kNMultiplicityBins = 15; | |
112 | Double_t multiplicityBins[kNMultiplicityBins+1]; | |
113 | for (Int_t i=0; i<=kNMultiplicityBins; i++) | |
114 | multiplicityBins[i] = -0.5 + i; | |
115 | multiplicityBins[kNMultiplicityBins] = 200; | |
bf58cbde | 116 | |
117 | // centrality | |
118 | const Int_t kNCentralityBins = 15; | |
119 | Double_t centralityBins[kNCentralityBins+1]; | |
120 | for (Int_t i=0; i<=kNCentralityBins; i++) | |
121 | centralityBins[i] = -0.5 + i * 500; | |
122 | ||
b1831bcb | 123 | // particle species |
124 | const Int_t kNSpeciesBins = 4; // pi, K, p, rest | |
125 | Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 }; | |
a75aacd6 | 126 | |
127 | trackBins[3] = multiplicityBins; | |
128 | iTrackBin[3] = kNMultiplicityBins; | |
129 | trackAxisTitle[3] = "multiplicity"; | |
130 | ||
131 | // selection depending on requested histogram | |
132 | Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead | |
133 | if (strcmp(reqHist, "NumberDensitypT") == 0) | |
134 | { | |
135 | axis = 0; | |
136 | title = "d^{2}N_{ch}/d#phid#eta"; | |
137 | } | |
138 | else if (strcmp(reqHist, "NumberDensityPhi") == 0) | |
139 | { | |
140 | axis = 1; | |
141 | title = "d^{2}N_{ch}/d#phid#eta"; | |
142 | } | |
bf58cbde | 143 | else if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0) |
144 | { | |
145 | axis = 2; | |
146 | title = "d^{2}N_{ch}/d#phid#eta"; | |
147 | } | |
a75aacd6 | 148 | else if (strcmp(reqHist, "SumpT") == 0) |
149 | { | |
150 | axis = 0; | |
151 | title = "d^{2}#Sigma p_{T}/d#phid#eta"; | |
152 | } | |
153 | else | |
154 | AliFatal(Form("Invalid histogram requested: %s", reqHist)); | |
155 | ||
156 | Int_t initRegions = fkRegions; | |
157 | ||
158 | if (axis == 0) | |
159 | { | |
160 | trackBins[2] = leadingpTBins; | |
161 | iTrackBin[2] = kNLeadingpTBins; | |
162 | trackAxisTitle[2] = "leading p_{T} (GeV/c)"; | |
163 | ||
164 | } | |
165 | else if (axis == 1) | |
166 | { | |
167 | nTrackVars = 5; | |
168 | initRegions = 1; | |
169 | ||
170 | iTrackBin[2] = kNLeadingpTBins2; | |
171 | trackBins[2] = leadingpTBins2; | |
172 | trackAxisTitle[2] = "leading p_{T} (GeV/c)"; | |
173 | ||
174 | iTrackBin[4] = kNLeadingPhiBins; | |
175 | trackBins[4] = leadingPhiBins; | |
bf58cbde | 176 | trackAxisTitle[4] = "#Delta #phi w.r.t. leading track"; |
177 | } | |
178 | else if (axis == 2) | |
179 | { | |
180 | nTrackVars = 5; | |
181 | initRegions = 1; | |
182 | ||
e0331fd9 | 183 | iTrackBin[0] = kNDeltaEtaBins; |
184 | trackBins[0] = deltaEtaBins; | |
185 | trackAxisTitle[0] = "#Delta#eta"; | |
186 | ||
bf58cbde | 187 | iTrackBin[2] = kNLeadingpTBins2; |
188 | trackBins[2] = leadingpTBins2; | |
189 | trackAxisTitle[2] = "leading p_{T} (GeV/c)"; | |
190 | ||
191 | trackBins[3] = centralityBins; | |
192 | iTrackBin[3] = kNCentralityBins; | |
193 | trackAxisTitle[3] = "centrality"; | |
194 | ||
195 | iTrackBin[4] = kNLeadingPhiBins; | |
196 | trackBins[4] = leadingPhiBins; | |
e0331fd9 | 197 | trackAxisTitle[4] = "#Delta#phi"; |
a75aacd6 | 198 | } |
199 | ||
200 | for (Int_t i=0; i<initRegions; i++) | |
201 | { | |
202 | fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin); | |
203 | ||
204 | for (Int_t j=0; j<nTrackVars; j++) | |
205 | { | |
206 | fTrackHist[i]->SetBinLimits(j, trackBins[j]); | |
207 | fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]); | |
208 | } | |
209 | ||
210 | SetStepNames(fTrackHist[i]); | |
211 | } | |
212 | ||
213 | // track 3rd and 4th axis --> event 1st and 2nd axis | |
214 | fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, 2, iTrackBin+2); | |
215 | ||
216 | fEventHist->SetBinLimits(0, trackBins[2]); | |
217 | fEventHist->SetVarTitle(0, trackAxisTitle[2]); | |
218 | ||
219 | fEventHist->SetBinLimits(1, trackBins[3]); | |
220 | fEventHist->SetVarTitle(1, trackAxisTitle[3]); | |
221 | ||
222 | SetStepNames(fEventHist); | |
b1831bcb | 223 | |
224 | iTrackBin[2] = kNSpeciesBins; | |
225 | ||
226 | fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 3, 3, iTrackBin); | |
227 | fTrackHistEfficiency->SetBinLimits(0, trackBins[0]); | |
228 | fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]); | |
229 | fTrackHistEfficiency->SetBinLimits(1, trackBins[1]); | |
230 | fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]); | |
231 | fTrackHistEfficiency->SetBinLimits(2, speciesBins); | |
232 | fTrackHistEfficiency->SetVarTitle(2, "particle species"); | |
a75aacd6 | 233 | } |
234 | ||
d1c75d06 | 235 | //_____________________________________________________________________________ |
236 | AliUEHist::AliUEHist(const AliUEHist &c) : | |
237 | TObject(), | |
238 | fkRegions(4), | |
239 | fEventHist(0), | |
b1831bcb | 240 | fTrackHistEfficiency(0), |
d1c75d06 | 241 | fEtaMin(0), |
242 | fEtaMax(0), | |
243 | fPtMin(0), | |
244 | fPtMax(0), | |
144bd037 | 245 | fContaminationEnhancement(0), |
d1c75d06 | 246 | fCombineMinMax(0), |
247 | fCache(0) | |
248 | { | |
249 | // | |
250 | // AliUEHist copy constructor | |
251 | // | |
252 | ||
253 | ((AliUEHist &) c).Copy(*this); | |
254 | } | |
255 | ||
a75aacd6 | 256 | //____________________________________________________________________ |
257 | void AliUEHist::SetStepNames(AliCFContainer* container) | |
258 | { | |
259 | // sets the names of the correction steps | |
260 | ||
261 | for (Int_t i=0; i<fgkCFSteps; i++) | |
262 | container->SetStepTitle(i, GetStepTitle((CFStep) i)); | |
263 | } | |
264 | ||
265 | //____________________________________________________________________ | |
266 | AliUEHist::~AliUEHist() | |
267 | { | |
268 | // Destructor | |
269 | ||
270 | for (Int_t i=0; i<fkRegions; i++) | |
271 | { | |
272 | if (fTrackHist[i]) | |
273 | { | |
274 | delete fTrackHist[i]; | |
275 | fTrackHist[i] = 0; | |
276 | } | |
277 | } | |
278 | ||
279 | if (fEventHist) | |
280 | { | |
281 | delete fEventHist; | |
282 | fEventHist = 0; | |
283 | } | |
284 | ||
b1831bcb | 285 | if (fTrackHistEfficiency) |
286 | { | |
287 | delete fTrackHistEfficiency; | |
288 | fTrackHistEfficiency = 0; | |
289 | } | |
290 | ||
a75aacd6 | 291 | if (fCache) |
292 | { | |
293 | delete fCache; | |
294 | fCache = 0; | |
295 | } | |
296 | } | |
297 | ||
298 | //____________________________________________________________________ | |
299 | AliUEHist &AliUEHist::operator=(const AliUEHist &c) | |
300 | { | |
301 | // assigment operator | |
302 | ||
303 | if (this != &c) | |
304 | ((AliUEHist &) c).Copy(*this); | |
305 | ||
306 | return *this; | |
307 | } | |
308 | ||
309 | //____________________________________________________________________ | |
310 | void AliUEHist::Copy(TObject& c) const | |
311 | { | |
312 | // copy function | |
313 | ||
314 | AliUEHist& target = (AliUEHist &) c; | |
315 | ||
316 | for (Int_t i=0; i<fkRegions; i++) | |
317 | if (fTrackHist[i]) | |
318 | target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone()); | |
319 | ||
320 | if (fEventHist) | |
321 | target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone()); | |
b1831bcb | 322 | |
323 | if (fTrackHistEfficiency) | |
324 | target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone()); | |
a75aacd6 | 325 | } |
326 | ||
327 | //____________________________________________________________________ | |
328 | Long64_t AliUEHist::Merge(TCollection* list) | |
329 | { | |
330 | // Merge a list of AliUEHist objects with this (needed for | |
331 | // PROOF). | |
332 | // Returns the number of merged objects (including this). | |
333 | ||
334 | if (!list) | |
335 | return 0; | |
336 | ||
337 | if (list->IsEmpty()) | |
338 | return 1; | |
339 | ||
340 | TIterator* iter = list->MakeIterator(); | |
341 | TObject* obj; | |
342 | ||
343 | // collections of objects | |
b1831bcb | 344 | const Int_t kMaxLists = fkRegions+2; |
a75aacd6 | 345 | TList** lists = new TList*[kMaxLists]; |
346 | ||
347 | for (Int_t i=0; i<kMaxLists; i++) | |
348 | lists[i] = new TList; | |
349 | ||
350 | Int_t count = 0; | |
351 | while ((obj = iter->Next())) { | |
352 | ||
353 | AliUEHist* entry = dynamic_cast<AliUEHist*> (obj); | |
354 | if (entry == 0) | |
355 | continue; | |
356 | ||
357 | for (Int_t i=0; i<fkRegions; i++) | |
358 | if (entry->fTrackHist[i]) | |
359 | lists[i]->Add(entry->fTrackHist[i]); | |
360 | ||
361 | lists[fkRegions]->Add(entry->fEventHist); | |
b1831bcb | 362 | lists[fkRegions+1]->Add(entry->fTrackHistEfficiency); |
a75aacd6 | 363 | |
364 | count++; | |
365 | } | |
366 | for (Int_t i=0; i<fkRegions; i++) | |
367 | if (fTrackHist[i]) | |
368 | fTrackHist[i]->Merge(lists[i]); | |
369 | ||
370 | fEventHist->Merge(lists[fkRegions]); | |
b1831bcb | 371 | fTrackHistEfficiency->Merge(lists[fkRegions+1]); |
a75aacd6 | 372 | |
373 | for (Int_t i=0; i<kMaxLists; i++) | |
374 | delete lists[i]; | |
375 | ||
376 | delete[] lists; | |
377 | ||
378 | return count+1; | |
379 | } | |
380 | ||
381 | //____________________________________________________________________ | |
382 | void AliUEHist::SetBinLimits(AliCFGridSparse* grid) | |
383 | { | |
384 | // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max | |
385 | ||
386 | if (fEtaMax > fEtaMin) | |
387 | grid->SetRangeUser(0, fEtaMin, fEtaMax); | |
388 | if (fPtMax > fPtMin) | |
389 | grid->SetRangeUser(1, fPtMin, fPtMax); | |
390 | } | |
391 | ||
392 | //____________________________________________________________________ | |
393 | void AliUEHist::ResetBinLimits(AliCFGridSparse* grid) | |
394 | { | |
395 | // resets all bin limits | |
396 | ||
397 | for (Int_t i=0; i<grid->GetNVar(); i++) | |
144bd037 | 398 | if (grid->GetGrid()->GetAxis(i)->TestBit(TAxis::kAxisRange)) |
a75aacd6 | 399 | grid->SetRangeUser(i, 0, -1); |
400 | } | |
401 | ||
144bd037 | 402 | //____________________________________________________________________ |
403 | void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax) | |
404 | { | |
405 | // prints the number of empty bins in the track end event histograms in the given step | |
406 | ||
407 | Int_t binBegin[4]; | |
408 | Int_t binEnd[4]; | |
409 | ||
410 | for (Int_t i=0; i<4; i++) | |
411 | { | |
412 | binBegin[i] = 1; | |
413 | binEnd[i] = fTrackHist[0]->GetNBins(i); | |
414 | } | |
415 | ||
416 | if (fEtaMax > fEtaMin) | |
417 | { | |
418 | binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin); | |
419 | binEnd[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax); | |
420 | } | |
421 | ||
422 | if (fPtMax > fPtMin) | |
423 | { | |
424 | binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin); | |
425 | binEnd[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax); | |
426 | } | |
427 | ||
428 | if (ptLeadMax > ptLeadMin) | |
429 | { | |
430 | binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin); | |
431 | binEnd[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax); | |
432 | } | |
433 | ||
434 | // start from multiplicity 1 | |
435 | binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1); | |
436 | ||
437 | for (Int_t region=0; region<fkRegions; region++) | |
438 | { | |
439 | Int_t total = 0; | |
440 | Int_t count = 0; | |
441 | Int_t vars[4]; | |
442 | ||
443 | for (Int_t i=0; i<4; i++) | |
444 | vars[i] = binBegin[i]; | |
445 | ||
446 | AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step); | |
447 | while (1) | |
448 | { | |
449 | if (grid->GetElement(vars) == 0) | |
450 | { | |
451 | Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f", | |
452 | grid->GetBinCenter(0, vars[0]), | |
453 | grid->GetBinCenter(1, vars[1]), | |
454 | grid->GetBinCenter(2, vars[2]), | |
455 | grid->GetBinCenter(3, vars[3]) | |
456 | ); | |
457 | count++; | |
458 | } | |
459 | ||
460 | vars[3]++; | |
461 | for (Int_t i=3; i>0; i--) | |
462 | { | |
463 | if (vars[i] == binEnd[i]+1) | |
464 | { | |
465 | vars[i] = binBegin[i]; | |
466 | vars[i-1]++; | |
467 | } | |
468 | } | |
469 | ||
470 | if (vars[0] == binEnd[0]+1) | |
471 | break; | |
472 | total++; | |
473 | } | |
474 | ||
475 | Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total); | |
476 | } | |
477 | } | |
478 | ||
a75aacd6 | 479 | //____________________________________________________________________ |
e0331fd9 | 480 | TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t twoD) |
a75aacd6 | 481 | { |
482 | // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events | |
483 | // | |
484 | // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead) | |
485 | // Histogram has to be deleted by the caller of the function | |
486 | ||
487 | // unzoom all axes | |
488 | ResetBinLimits(fTrackHist[region]->GetGrid(step)); | |
489 | ResetBinLimits(fEventHist->GetGrid(step)); | |
490 | ||
491 | SetBinLimits(fTrackHist[region]->GetGrid(step)); | |
492 | ||
493 | TH1D* tracks = 0; | |
494 | ||
495 | if (ptLeadMin < 0) | |
496 | { | |
497 | tracks = fTrackHist[region]->ShowProjection(2, step); | |
498 | tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle()); | |
499 | if (fCombineMinMax && region == kMin) | |
500 | { | |
501 | ResetBinLimits(fTrackHist[kMax]->GetGrid(step)); | |
502 | SetBinLimits(fTrackHist[kMax]->GetGrid(step)); | |
503 | ||
504 | TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step); | |
505 | tracks->Add(tracks2); | |
506 | ||
507 | ResetBinLimits(fTrackHist[kMax]->GetGrid(step)); | |
508 | } | |
509 | ||
510 | // normalize to get a density (deta dphi) | |
511 | TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0); | |
512 | Float_t phiRegion = TMath::TwoPi() / 3; | |
513 | if (!fCombineMinMax && region == kMin) | |
514 | phiRegion /= 2; | |
515 | Float_t normalization = phiRegion * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst())); | |
516 | //Printf("Normalization: %f", normalization); | |
517 | tracks->Scale(1.0 / normalization); | |
518 | ||
519 | TH1D* events = fEventHist->ShowProjection(0, step); | |
520 | tracks->Divide(events); | |
521 | } | |
522 | else | |
523 | { | |
e741fadd | 524 | // 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 |
525 | // therefore the number density must be calculated here per leading pT bin and then summed | |
526 | ||
e0331fd9 | 527 | if (multBinEnd > multBinBegin) |
528 | Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd); | |
bf58cbde | 529 | fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd); |
530 | fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd); | |
531 | ||
e741fadd | 532 | Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin); |
533 | Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax); | |
a75aacd6 | 534 | |
bf58cbde | 535 | TH1D* events = fEventHist->ShowProjection(0, step); |
536 | ||
e741fadd | 537 | for (Int_t bin=firstBin; bin<=lastBin; bin++) |
538 | { | |
539 | fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin); | |
e0331fd9 | 540 | TH1D* tracksTmp = 0; |
541 | if (twoD) | |
542 | tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0); | |
543 | else | |
544 | tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4); | |
e741fadd | 545 | Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral()); |
546 | fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1); | |
547 | ||
548 | // normalize to get a density (deta dphi) | |
549 | TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0); | |
550 | Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1) * (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst())); | |
551 | //Printf("Normalization: %f", normalization); | |
552 | tracksTmp->Scale(1.0 / normalization); | |
553 | ||
e741fadd | 554 | Int_t nEvents = (Int_t) events->GetBinContent(bin); |
555 | if (nEvents > 0) | |
556 | tracksTmp->Scale(1.0 / nEvents); | |
557 | Printf("Calculated histogram in bin %d --> %d events", bin, nEvents); | |
558 | ||
559 | if (!tracks) | |
560 | tracks = tracksTmp; | |
561 | else | |
562 | { | |
563 | tracks->Add(tracksTmp); | |
564 | delete tracksTmp; | |
565 | } | |
566 | } | |
bf58cbde | 567 | |
568 | delete events; | |
a75aacd6 | 569 | } |
570 | ||
571 | ResetBinLimits(fTrackHist[region]->GetGrid(step)); | |
bf58cbde | 572 | ResetBinLimits(fEventHist->GetGrid(step)); |
a75aacd6 | 573 | |
574 | return tracks; | |
575 | } | |
576 | ||
577 | //____________________________________________________________________ | |
578 | void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2) | |
579 | { | |
580 | // corrects from step1 to step2 by multiplying the tracks with trackCorrection | |
581 | // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4) | |
582 | // if var2 >= 0 a two dimension correction is assumed in trackCorrection | |
583 | // | |
584 | // if trackCorrection is 0, just copies content from step1 to step2 | |
585 | ||
586 | for (Int_t region=0; region<fkRegions; region++) | |
144bd037 | 587 | CorrectTracks(step1, step2, region, trackCorrection, var1, var2); |
588 | } | |
589 | ||
590 | //____________________________________________________________________ | |
591 | void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2) | |
592 | { | |
593 | // | |
594 | // see documentation of CorrectTracks above | |
595 | // | |
596 | ||
597 | if (!fTrackHist[region]) | |
598 | return; | |
5a53c6f2 | 599 | |
144bd037 | 600 | THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid(); |
601 | THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid(); | |
602 | ||
603 | // clear target histogram | |
604 | target->Reset(); | |
605 | ||
606 | if (trackCorrection != 0) | |
a75aacd6 | 607 | { |
144bd037 | 608 | if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX()) |
609 | AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX())); | |
a75aacd6 | 610 | |
144bd037 | 611 | if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY()) |
612 | AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY())); | |
613 | } | |
614 | ||
615 | // optimized implementation | |
616 | for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++) | |
617 | { | |
618 | Int_t bins[5]; | |
619 | Double_t value = grid->GetBinContent(binIdx, bins); | |
620 | Double_t error = grid->GetBinError(binIdx); | |
a75aacd6 | 621 | |
622 | if (trackCorrection != 0) | |
623 | { | |
144bd037 | 624 | if (var2 < 0) |
a75aacd6 | 625 | { |
144bd037 | 626 | value *= trackCorrection->GetBinContent(bins[var1]); |
627 | error *= trackCorrection->GetBinContent(bins[var1]); | |
628 | } | |
629 | else | |
630 | { | |
631 | value *= trackCorrection->GetBinContent(bins[var1], bins[var2]); | |
632 | error *= trackCorrection->GetBinContent(bins[var1], bins[var2]); | |
a75aacd6 | 633 | } |
a75aacd6 | 634 | } |
144bd037 | 635 | |
636 | target->SetBinContent(bins, value); | |
637 | target->SetBinError(bins, error); | |
a75aacd6 | 638 | } |
5a53c6f2 | 639 | |
640 | Printf("AliUEHist::CorrectTracks: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0); | |
a75aacd6 | 641 | } |
642 | ||
643 | //____________________________________________________________________ | |
644 | void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1D* eventCorrection, Int_t var) | |
645 | { | |
646 | // corrects from step1 to step2 by multiplying the events with eventCorrection | |
647 | // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1) | |
648 | // | |
649 | // if eventCorrection is 0, just copies content from step1 to step2 | |
650 | ||
651 | AliCFGridSparse* grid = fEventHist->GetGrid(step1); | |
652 | AliCFGridSparse* target = fEventHist->GetGrid(step2); | |
653 | ||
144bd037 | 654 | // clear target histogram |
655 | target->GetGrid()->Reset(); | |
656 | ||
a75aacd6 | 657 | if (eventCorrection != 0 && grid->GetNBins(var) != eventCorrection->GetNbinsX()) |
658 | AliFatal(Form("Invalid binning: %d %d", grid->GetNBins(var), eventCorrection->GetNbinsX())); | |
659 | ||
660 | Int_t bins[2]; | |
661 | for (Int_t x = 1; x <= grid->GetNBins(0); x++) | |
662 | { | |
663 | bins[0] = x; | |
664 | for (Int_t y = 1; y <= grid->GetNBins(1); y++) | |
665 | { | |
666 | bins[1] = y; | |
667 | ||
668 | Double_t value = grid->GetElement(bins); | |
669 | if (value != 0) | |
670 | { | |
671 | Double_t error = grid->GetElementError(bins); | |
672 | ||
673 | if (eventCorrection != 0) | |
674 | { | |
675 | value *= eventCorrection->GetBinContent(bins[var]); | |
676 | error *= eventCorrection->GetBinContent(bins[var]); | |
677 | } | |
678 | ||
679 | target->SetElement(bins, value); | |
680 | target->SetElementError(bins, error); | |
681 | } | |
682 | } | |
683 | } | |
5a53c6f2 | 684 | |
685 | Printf("AliUEHist::CorrectEvents: Corrected from %f to %f entries. Correction histogram: %f entries (integral: %f)", grid->GetEntries(), target->GetEntries(), (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0); | |
a75aacd6 | 686 | } |
687 | ||
688 | //____________________________________________________________________ | |
689 | void AliUEHist::Correct(AliUEHist* corrections) | |
690 | { | |
691 | // applies the given corrections to extract from the step kCFStepReconstructed all previous steps | |
692 | // | |
693 | // in this object the data is expected in the step kCFStepReconstructed | |
694 | ||
695 | // ---- track level | |
696 | ||
697 | // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used) | |
698 | // extracted as function of leading pT | |
144bd037 | 699 | for (Int_t region = 0; region < fkRegions; region++) |
700 | { | |
701 | if (!fTrackHist[region]) | |
702 | continue; | |
ada1a03f | 703 | |
704 | const char* projAxis = "z"; | |
705 | Int_t secondBin = -1; | |
706 | ||
707 | if (fTrackHist[region]->GetNVar() == 5) | |
708 | { | |
709 | projAxis = "yz"; | |
710 | secondBin = 4; | |
711 | } | |
712 | ||
713 | #if 0 | |
714 | TH1* leadingBias = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis); // from MC | |
715 | Printf("WARNING: Using MC bias correction"); | |
716 | #else | |
717 | TH1* leadingBias = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis); // from data | |
718 | #endif | |
719 | CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBias, 2, secondBin); | |
144bd037 | 720 | if (region == kMin && fCombineMinMax) |
721 | { | |
ada1a03f | 722 | CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBias, 2, secondBin); |
144bd037 | 723 | delete leadingBias; |
724 | break; | |
725 | } | |
726 | delete leadingBias; | |
727 | } | |
a75aacd6 | 728 | CorrectEvents(kCFStepReconstructed, kCFStepTracked, 0, 0); |
a75aacd6 | 729 | |
730 | // correct with kCFStepTracked --> kCFStepTrackedOnlyPrim | |
731 | TH2D* contamination = corrections->GetTrackingContamination(); | |
144bd037 | 732 | if (corrections->fContaminationEnhancement) |
733 | { | |
734 | Printf("Applying contamination enhancement"); | |
735 | ||
736 | for (Int_t x=1; x<=contamination->GetNbinsX(); x++) | |
737 | for (Int_t y=1; y<=contamination->GetNbinsX(); y++) | |
738 | contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y)))); | |
739 | } | |
a75aacd6 | 740 | CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1); |
741 | CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0); | |
742 | delete contamination; | |
743 | ||
144bd037 | 744 | // correct with kCFStepTrackedOnlyPrim --> kCFStepAnaTopology |
745 | TH2D* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection(); | |
746 | CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0, 1); | |
747 | CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 0); | |
748 | delete efficiencyCorrection; | |
a75aacd6 | 749 | |
750 | // copy | |
751 | CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1); | |
752 | CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0); | |
753 | ||
754 | // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor | |
755 | // practically independent of low pT cut | |
756 | TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1); | |
757 | ||
758 | // convert stage from true multiplicity to observed multiplicity by simple conversion factor | |
759 | TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2"); | |
760 | vertexCorrectionObs->Reset(); | |
761 | ||
144bd037 | 762 | TF1* func = new TF1("func", "[1]+[0]/(x-[2])"); |
9217ab56 | 763 | // some defaults |
764 | func->SetParameters(0.1, 1, -0.7); | |
e741fadd | 765 | vertexCorrection->Fit(func, "0I", "", 0, 3); |
a75aacd6 | 766 | for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++) |
144bd037 | 767 | { |
768 | Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i); | |
e741fadd | 769 | if (xPos < 3) |
144bd037 | 770 | vertexCorrectionObs->SetBinContent(i, func->Eval(xPos)); |
771 | else | |
772 | vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos)); | |
773 | } | |
a75aacd6 | 774 | |
9217ab56 | 775 | #if 1 |
144bd037 | 776 | new TCanvas; |
a75aacd6 | 777 | vertexCorrection->DrawCopy(); |
778 | vertexCorrectionObs->SetLineColor(2); | |
144bd037 | 779 | vertexCorrectionObs->DrawCopy("same"); |
ada1a03f | 780 | func->SetRange(0, 4); |
781 | func->DrawClone("same"); | |
b1831bcb | 782 | #endif |
a75aacd6 | 783 | |
784 | CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3); | |
785 | CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1); | |
786 | delete vertexCorrectionObs; | |
787 | delete vertexCorrection; | |
ada1a03f | 788 | delete func; |
a75aacd6 | 789 | |
790 | // copy | |
791 | CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1); | |
792 | CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0); | |
793 | } | |
794 | ||
795 | //____________________________________________________________________ | |
b1831bcb | 796 | TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source) |
a75aacd6 | 797 | { |
798 | // creates a track-level efficiency by dividing step2 by step1 | |
799 | // projected to axis1 and axis2 (optional if >= 0) | |
b1831bcb | 800 | // |
801 | // source: 0 = fTrackHist; 1 = fTrackHistEfficiency | |
a75aacd6 | 802 | |
803 | // integrate over regions | |
804 | // cache it for efficiency (usually more than one efficiency is requested) | |
b1831bcb | 805 | |
806 | AliCFContainer* sourceContainer = 0; | |
807 | ||
808 | if (source == 0) | |
a75aacd6 | 809 | { |
b1831bcb | 810 | if (!fCache) |
811 | { | |
812 | fCache = (AliCFContainer*) fTrackHist[0]->Clone(); | |
813 | for (Int_t i = 1; i < fkRegions; i++) | |
814 | if (fTrackHist[i]) | |
815 | fCache->Add(fTrackHist[i]); | |
816 | } | |
817 | sourceContainer = fCache; | |
a75aacd6 | 818 | } |
b1831bcb | 819 | else if (source == 1) |
820 | { | |
821 | sourceContainer = fTrackHistEfficiency; | |
822 | // step offset because we start with kCFStepAnaTopology | |
823 | step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology); | |
824 | step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology); | |
825 | } | |
826 | else | |
827 | return 0; | |
828 | ||
a75aacd6 | 829 | // reset all limits and set the right ones except those in axis1 and axis2 |
b1831bcb | 830 | ResetBinLimits(sourceContainer->GetGrid(step1)); |
831 | ResetBinLimits(sourceContainer->GetGrid(step2)); | |
a75aacd6 | 832 | if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0) |
833 | { | |
b1831bcb | 834 | sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax); |
835 | sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax); | |
a75aacd6 | 836 | } |
837 | if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1) | |
838 | { | |
b1831bcb | 839 | sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax); |
840 | sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax); | |
a75aacd6 | 841 | } |
842 | ||
843 | TH1* measured = 0; | |
844 | TH1* generated = 0; | |
845 | ||
144bd037 | 846 | if (axis2 >= 0) |
a75aacd6 | 847 | { |
b1831bcb | 848 | generated = sourceContainer->Project(axis1, axis2, step1); |
849 | measured = sourceContainer->Project(axis1, axis2, step2); | |
a75aacd6 | 850 | } |
851 | else | |
852 | { | |
b1831bcb | 853 | generated = sourceContainer->Project(axis1, step1); |
854 | measured = sourceContainer->Project(axis1, step2); | |
a75aacd6 | 855 | } |
856 | ||
144bd037 | 857 | // check for bins with less than 100 entries, print warning |
858 | Int_t binBegin[2]; | |
859 | Int_t binEnd[2]; | |
860 | ||
861 | binBegin[0] = 1; | |
862 | binBegin[1] = 1; | |
863 | ||
864 | binEnd[0] = generated->GetNbinsX(); | |
865 | binEnd[1] = generated->GetNbinsY(); | |
866 | ||
867 | if (fEtaMax > fEtaMin) | |
868 | { | |
869 | if (axis1 == 0) | |
870 | { | |
871 | binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin); | |
872 | binEnd[0] = generated->GetXaxis()->FindBin(fEtaMax); | |
873 | } | |
874 | if (axis2 == 0) | |
875 | { | |
876 | binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin); | |
877 | binEnd[1] = generated->GetYaxis()->FindBin(fEtaMax); | |
878 | } | |
879 | } | |
880 | ||
881 | if (fPtMax > fPtMin) | |
882 | { | |
883 | // TODO this is just checking up to 15 for now | |
884 | Float_t ptMax = TMath::Min((Float_t) 15., fPtMax); | |
885 | if (axis1 == 1) | |
886 | { | |
887 | binBegin[0] = generated->GetXaxis()->FindBin(fPtMin); | |
888 | binEnd[0] = generated->GetXaxis()->FindBin(ptMax); | |
889 | } | |
890 | if (axis2 == 1) | |
891 | { | |
892 | binBegin[1] = generated->GetYaxis()->FindBin(fPtMin); | |
893 | binEnd[1] = generated->GetYaxis()->FindBin(ptMax); | |
894 | } | |
895 | } | |
896 | ||
897 | Int_t total = 0; | |
898 | Int_t count = 0; | |
899 | Int_t vars[2]; | |
900 | ||
901 | for (Int_t i=0; i<2; i++) | |
902 | vars[i] = binBegin[i]; | |
903 | ||
904 | const Int_t limit = 50; | |
905 | while (1) | |
906 | { | |
907 | if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit) | |
908 | { | |
909 | Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0])); | |
910 | count++; | |
911 | } | |
912 | else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit) | |
913 | { | |
914 | Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)", | |
915 | generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), | |
916 | generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]), | |
917 | generated->GetBinContent(vars[0], vars[1])); | |
918 | count++; | |
919 | } | |
920 | ||
921 | vars[1]++; | |
922 | if (vars[1] == binEnd[1]+1) | |
923 | { | |
924 | vars[1] = binBegin[1]; | |
925 | vars[0]++; | |
926 | } | |
927 | ||
928 | if (vars[0] == binEnd[0]+1) | |
929 | break; | |
930 | total++; | |
931 | } | |
932 | ||
933 | Printf("Correction has %d empty bins (out of %d bins)", count, total); | |
934 | ||
935 | measured->Divide(measured, generated, 1, 1, "B"); | |
a75aacd6 | 936 | |
937 | delete generated; | |
938 | ||
b1831bcb | 939 | ResetBinLimits(sourceContainer->GetGrid(step1)); |
940 | ResetBinLimits(sourceContainer->GetGrid(step2)); | |
941 | ||
a75aacd6 | 942 | return measured; |
943 | } | |
944 | ||
945 | //____________________________________________________________________ | |
946 | TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax) | |
947 | { | |
948 | // creates a event-level efficiency by dividing step2 by step1 | |
949 | // projected to axis1 and axis2 (optional if >= 0) | |
950 | ||
951 | if (ptLeadMax > ptLeadMin) | |
952 | { | |
953 | fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax); | |
954 | fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax); | |
955 | } | |
956 | ||
957 | TH1* measured = 0; | |
958 | TH1* generated = 0; | |
959 | ||
144bd037 | 960 | if (axis2 >= 0) |
a75aacd6 | 961 | { |
962 | generated = fEventHist->Project(axis1, axis2, step1); | |
963 | measured = fEventHist->Project(axis1, axis2, step2); | |
964 | } | |
965 | else | |
966 | { | |
967 | generated = fEventHist->Project(axis1, step1); | |
968 | measured = fEventHist->Project(axis1, step2); | |
969 | } | |
970 | ||
144bd037 | 971 | measured->Divide(measured, generated, 1, 1, "B"); |
a75aacd6 | 972 | |
973 | delete generated; | |
974 | ||
975 | if (ptLeadMax > ptLeadMin) | |
976 | { | |
977 | fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1); | |
978 | fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1); | |
979 | } | |
980 | ||
981 | return measured; | |
982 | } | |
983 | ||
984 | //____________________________________________________________________ | |
985 | void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2) | |
986 | { | |
987 | // weights each entry of the 3d histogram hist1 with the 1d histogram hist2 | |
988 | // where the matching is done of the z axis of hist1 with the x axis of hist2 | |
989 | ||
990 | if (hist1->GetNbinsZ() != hist2->GetNbinsX()) | |
991 | AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX())); | |
992 | ||
993 | for (Int_t x=1; x<=hist1->GetNbinsX(); x++) | |
994 | { | |
995 | for (Int_t y=1; y<=hist1->GetNbinsY(); y++) | |
996 | { | |
997 | for (Int_t z=1; z<=hist1->GetNbinsZ(); z++) | |
998 | { | |
999 | if (hist2->GetBinContent(z) > 0) | |
1000 | { | |
1001 | hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z)); | |
1002 | hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z)); | |
1003 | } | |
1004 | else | |
1005 | { | |
1006 | hist1->SetBinContent(x, y, z, 0); | |
1007 | hist1->SetBinError(x, y, z, 0); | |
1008 | } | |
1009 | } | |
1010 | } | |
1011 | } | |
1012 | } | |
1013 | ||
1014 | //____________________________________________________________________ | |
144bd037 | 1015 | TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax) |
a75aacd6 | 1016 | { |
1017 | // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1) | |
144bd037 | 1018 | // in the given region (sum over all regions is calculated if region == -1) |
a75aacd6 | 1019 | // done by weighting the track-level distribution with the number of events as function of leading pT |
1020 | // and then calculating the ratio between the distributions | |
1021 | // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx" | |
1022 | // no projection is done if axis == 0 | |
1023 | ||
144bd037 | 1024 | AliCFContainer* tmp = 0; |
1025 | ||
1026 | if (region == -1) | |
1027 | { | |
1028 | tmp = (AliCFContainer*) fTrackHist[0]->Clone(); | |
1029 | for (Int_t i = 1; i < fkRegions; i++) | |
1030 | if (fTrackHist[i]) | |
1031 | tmp->Add(fTrackHist[i]); | |
1032 | } | |
1033 | else if (region == kMin && fCombineMinMax) | |
1034 | { | |
1035 | tmp = (AliCFContainer*) fTrackHist[kMin]->Clone(); | |
1036 | tmp->Add(fTrackHist[kMax]); | |
1037 | } | |
1038 | else | |
1039 | tmp = fTrackHist[region]; | |
1040 | ||
1041 | ResetBinLimits(tmp->GetGrid(step1)); | |
1042 | ResetBinLimits(fEventHist->GetGrid(step1)); | |
1043 | SetBinLimits(tmp->GetGrid(step1)); | |
1044 | ||
1045 | ResetBinLimits(tmp->GetGrid(step2)); | |
1046 | ResetBinLimits(fEventHist->GetGrid(step2)); | |
1047 | SetBinLimits(tmp->GetGrid(step2)); | |
a75aacd6 | 1048 | |
df269636 | 1049 | TH1D* events1 = (TH1D*)fEventHist->Project(0, step1); |
1050 | TH3D* hist1 = (TH3D*)tmp->Project(0, tmp->GetNVar()-1, 2, step1); | |
a75aacd6 | 1051 | WeightHistogram(hist1, events1); |
1052 | ||
df269636 | 1053 | TH1D* events2 = (TH1D*)fEventHist->Project(0, step2); |
1054 | TH3D* hist2 = (TH3D*)tmp->Project(0, tmp->GetNVar()-1, 2, step2); | |
a75aacd6 | 1055 | WeightHistogram(hist2, events2); |
1056 | ||
1057 | TH1* generated = hist1; | |
1058 | TH1* measured = hist2; | |
1059 | ||
1060 | if (axis) | |
1061 | { | |
1062 | if (leadPtMax > leadPtMin) | |
1063 | { | |
1064 | hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax); | |
1065 | hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax); | |
1066 | } | |
1067 | ||
1068 | if (fEtaMax > fEtaMin && !TString(axis).Contains("x")) | |
1069 | { | |
1070 | hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax); | |
1071 | hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax); | |
1072 | } | |
ada1a03f | 1073 | |
a75aacd6 | 1074 | generated = hist1->Project3D(axis); |
1075 | measured = hist2->Project3D(axis); | |
1076 | ||
1077 | // delete hists here if projection has been done | |
1078 | delete hist1; | |
1079 | delete hist2; | |
1080 | } | |
1081 | ||
1082 | measured->Divide(generated); | |
1083 | ||
1084 | delete events1; | |
1085 | delete events2; | |
1086 | delete generated; | |
144bd037 | 1087 | |
1088 | ResetBinLimits(tmp->GetGrid(step1)); | |
1089 | ResetBinLimits(tmp->GetGrid(step2)); | |
1090 | ||
29e8465c | 1091 | if ((region == -1) || (region == kMin && fCombineMinMax)) |
144bd037 | 1092 | delete tmp; |
a75aacd6 | 1093 | |
1094 | return measured; | |
1095 | } | |
1096 | ||
1097 | //____________________________________________________________________ | |
1098 | TH2D* AliUEHist::GetTrackingEfficiency() | |
1099 | { | |
1100 | // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim | |
1101 | // integrates over the regions and all other variables than pT and eta to increase the statistics | |
1102 | // | |
1103 | // returned histogram has to be deleted by the user | |
1104 | ||
1105 | return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1)); | |
1106 | } | |
1107 | ||
1108 | //____________________________________________________________________ | |
1109 | TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis) | |
1110 | { | |
1111 | // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim | |
1112 | // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics | |
1113 | ||
1114 | return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis)); | |
1115 | } | |
1116 | ||
1117 | //____________________________________________________________________ | |
1118 | TH2D* AliUEHist::GetTrackingCorrection() | |
1119 | { | |
1120 | // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked | |
1121 | // integrates over the regions and all other variables than pT and eta to increase the statistics | |
1122 | // | |
1123 | // returned histogram has to be deleted by the user | |
1124 | ||
1125 | return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1)); | |
1126 | } | |
1127 | ||
1128 | //____________________________________________________________________ | |
1129 | TH1D* AliUEHist::GetTrackingCorrection(Int_t axis) | |
1130 | { | |
1131 | // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked | |
1132 | // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics | |
1133 | ||
1134 | return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis)); | |
1135 | } | |
1136 | ||
144bd037 | 1137 | //____________________________________________________________________ |
1138 | TH2D* AliUEHist::GetTrackingEfficiencyCorrection() | |
1139 | { | |
1140 | // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked | |
1141 | // integrates over the regions and all other variables than pT and eta to increase the statistics | |
1142 | // | |
1143 | // returned histogram has to be deleted by the user | |
1144 | ||
1145 | return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1)); | |
1146 | } | |
1147 | ||
1148 | //____________________________________________________________________ | |
1149 | TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis) | |
1150 | { | |
1151 | // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked | |
1152 | // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics | |
1153 | ||
1154 | return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis)); | |
1155 | } | |
1156 | ||
a75aacd6 | 1157 | //____________________________________________________________________ |
1158 | TH2D* AliUEHist::GetTrackingContamination() | |
1159 | { | |
1160 | // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked | |
1161 | // integrates over the regions and all other variables than pT and eta to increase the statistics | |
1162 | // | |
1163 | // returned histogram has to be deleted by the user | |
1164 | ||
1165 | return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1)); | |
1166 | } | |
1167 | ||
1168 | //____________________________________________________________________ | |
1169 | TH1D* AliUEHist::GetTrackingContamination(Int_t axis) | |
1170 | { | |
1171 | // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked | |
1172 | // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics | |
1173 | ||
1174 | return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis)); | |
1175 | } | |
1176 | ||
1177 | //____________________________________________________________________ | |
1178 | const char* AliUEHist::GetRegionTitle(Region region) | |
1179 | { | |
1180 | // returns the name of the given region | |
1181 | ||
1182 | switch (region) | |
1183 | { | |
1184 | case kToward: | |
1185 | return "Towards"; | |
1186 | case kAway: | |
1187 | return "Away"; | |
1188 | case kMin: | |
1189 | return (fCombineMinMax) ? "Transverse" : "Min"; | |
1190 | case kMax: | |
1191 | return "Max"; | |
1192 | } | |
1193 | ||
1194 | return 0; | |
1195 | } | |
1196 | ||
1197 | //____________________________________________________________________ | |
1198 | const char* AliUEHist::GetStepTitle(CFStep step) | |
1199 | { | |
1200 | // returns the name of the given step | |
1201 | ||
1202 | switch (step) | |
1203 | { | |
1204 | case kCFStepAll: | |
1205 | return "All events"; | |
1206 | case kCFStepTriggered: | |
1207 | return "Triggered"; | |
1208 | case kCFStepVertex: | |
1209 | return "Primary Vertex"; | |
1210 | case kCFStepAnaTopology: | |
1211 | return "Required analysis topology"; | |
1212 | case kCFStepTrackedOnlyPrim: | |
1213 | return "Tracked (matched MC, only primaries)"; | |
1214 | case kCFStepTracked: | |
1215 | return "Tracked (matched MC, all)"; | |
1216 | case kCFStepReconstructed: | |
1217 | return "Reconstructed"; | |
1218 | case kCFStepRealLeading: | |
1219 | return "Correct leading particle identified"; | |
1220 | case kCFStepBiasStudy: | |
1221 | return "Bias study applying tracking efficiency"; | |
1222 | case kCFStepBiasStudy2: | |
1223 | return "Bias study applying tracking efficiency in two steps"; | |
1224 | } | |
1225 | ||
1226 | return 0; | |
1227 | } | |
b1831bcb | 1228 | |
1229 | //____________________________________________________________________ | |
1230 | void AliUEHist::CopyReconstructedData(AliUEHist* from) | |
1231 | { | |
1232 | // copies those histograms extracted from ESD to this object | |
1233 | ||
1234 | // TODO at present only the pointers are copied | |
1235 | ||
1236 | for (Int_t region=0; region<4; region++) | |
1237 | { | |
1238 | if (!fTrackHist[region]) | |
1239 | continue; | |
1240 | ||
1241 | fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed)); | |
df269636 | 1242 | //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim)); |
b1831bcb | 1243 | fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy)); |
1244 | } | |
1245 | ||
1246 | fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed)); | |
df269636 | 1247 | //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim)); |
b1831bcb | 1248 | fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy, from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy)); |
1249 | } | |
6f803f6c | 1250 | |
1251 | //____________________________________________________________________ | |
df269636 | 1252 | void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose) |
6f803f6c | 1253 | { |
1254 | // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency | |
1255 | ||
1256 | Float_t fitRangeBegin = 5.01; | |
1257 | Float_t fitRangeEnd = 14.99; | |
1258 | Float_t extendRangeBegin = 10.01; | |
1259 | ||
1260 | TH1* obj = GetTrackingEfficiency(1); | |
1261 | ||
df269636 | 1262 | if (verbose) |
1263 | new TCanvas; obj->Draw(); | |
1264 | obj->Fit("pol0", (verbose) ? "+" : "+0", "SAME", fitRangeBegin, fitRangeEnd); | |
6f803f6c | 1265 | |
1266 | Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0); | |
1267 | ||
1268 | obj = GetTrackingContamination(1); | |
1269 | ||
df269636 | 1270 | if (verbose) |
1271 | new TCanvas; obj->Draw(); | |
1272 | obj->Fit("pol0", (verbose) ? "+" : "+0", "SAME", fitRangeBegin, fitRangeEnd); | |
6f803f6c | 1273 | |
1274 | Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0); | |
1275 | ||
1276 | 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); | |
1277 | ||
1278 | // extend up to pT 100 | |
1279 | for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++) | |
1280 | for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++) | |
df269636 | 1281 | for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis |
6f803f6c | 1282 | { |
6f803f6c | 1283 | Int_t bins[3]; |
1284 | bins[0] = x; | |
1285 | bins[1] = y; | |
1286 | bins[2] = z; | |
1287 | ||
df269636 | 1288 | fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100); |
1289 | fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff); | |
1290 | fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont); | |
6f803f6c | 1291 | } |
1292 | } | |
e741fadd | 1293 | |
1294 | void AliUEHist::AdditionalDPhiCorrection(Int_t step) | |
1295 | { | |
1296 | // corrects the dphi distribution with an extra factor close to dphi ~ 0 | |
1297 | ||
1298 | Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection."); | |
1299 | ||
1300 | THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid(); | |
1301 | ||
1302 | // optimized implementation | |
1303 | for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++) | |
1304 | { | |
1305 | Int_t bins[5]; | |
1306 | Double_t value = grid->GetBinContent(binIdx, bins); | |
1307 | Double_t error = grid->GetBinError(binIdx); | |
1308 | ||
1309 | Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]); | |
1310 | if (TMath::Abs(binCenter) < 0.2) | |
1311 | { | |
1312 | value *= 0.985; | |
1313 | error *= 0.985; | |
1314 | } | |
1315 | else if (TMath::Abs(binCenter) < 0.3) | |
1316 | { | |
1317 | value *= 0.9925; | |
1318 | error *= 0.9925; | |
1319 | } | |
1320 | ||
1321 | grid->SetBinContent(bins, value); | |
1322 | grid->SetBinError(bins, error); | |
1323 | } | |
1324 | } |