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