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