]>
Commit | Line | Data |
---|---|---|
1 | #include "AliPoissonCalculator.h" | |
2 | #include "AliForwardCorrectionManager.h" | |
3 | #include <TH2D.h> | |
4 | #include <TBrowser.h> | |
5 | #include <TROOT.h> | |
6 | #include <TMath.h> | |
7 | #include <TList.h> | |
8 | #include <iostream> | |
9 | #include <TAxis.h> | |
10 | ||
11 | // | |
12 | // A class to calculate the multiplicity in @f$(x,y)@f$ bins | |
13 | // using Poisson statistics. | |
14 | // | |
15 | // The input is assumed to be binned in @f$(x,y)@f$ as described by | |
16 | // the 2D histogram passwd to the Reset member function. | |
17 | // | |
18 | // The data is grouped in to regions as defined by the parameters | |
19 | // fXLumping and fYLumping. The total number of cells and number of | |
20 | // empty cells is then calculate in each region. The mean | |
21 | // multiplicity over the region is then determined as | |
22 | // | |
23 | // @f[ | |
24 | // \lange m\rangle = -\log\left(\frac{e}{t}\right) | |
25 | // @f] | |
26 | // where @f$ e@f$ is the number of empty cells and @f$t@f$ is the | |
27 | // total number of cells in the region. A correction for counting | |
28 | // statistics, is then applied | |
29 | // @f{eqnarray*} | |
30 | // c &=& \frac{1}{1 - \exp{-\lange m\rangle}} | |
31 | // &=& \frac{1}{1 - \frac{e}{t}} | |
32 | // @f{eqnarray*} | |
33 | // and the final number in each cell is then | |
34 | // @f[ | |
35 | // h_i c \lange m\rangle | |
36 | // @f] | |
37 | // where @f$h_i@f$ is the number of hits in the cell @f$i@f$ | |
38 | // | |
39 | // | |
40 | ||
41 | namespace { | |
42 | const char* kBasicN = "basic"; | |
43 | const char* kEmptyN = "empty"; | |
44 | const char* kTotalN = "total"; | |
45 | const char* kBasicT = "Basic number of hits"; | |
46 | const char* kEmptyT = "Empty number of bins/region"; | |
47 | const char* kTotalT = "Total number of bins/region"; | |
48 | } | |
49 | ||
50 | ||
51 | ||
52 | ||
53 | //____________________________________________________________________ | |
54 | AliPoissonCalculator::AliPoissonCalculator() | |
55 | : TNamed(), | |
56 | fXLumping(32), | |
57 | fYLumping(4), | |
58 | fTotal(0), | |
59 | fEmpty(0), | |
60 | fBasic(0), | |
61 | fEmptyVsTotal(0), | |
62 | fMean(0), | |
63 | fOcc(0), | |
64 | fCorr(0) | |
65 | { | |
66 | // | |
67 | // CTOR | |
68 | // | |
69 | } | |
70 | ||
71 | //____________________________________________________________________ | |
72 | AliPoissonCalculator::AliPoissonCalculator(const char*) | |
73 | : TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"), | |
74 | fXLumping(32), | |
75 | fYLumping(4), | |
76 | fTotal(0), | |
77 | fEmpty(0), | |
78 | fBasic(0), | |
79 | fEmptyVsTotal(0), | |
80 | fMean(0), | |
81 | fOcc(0), | |
82 | fCorr(0) | |
83 | { | |
84 | // | |
85 | // CTOR | |
86 | // | |
87 | } | |
88 | //____________________________________________________________________ | |
89 | AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o) | |
90 | : TNamed(o), | |
91 | fXLumping(o.fXLumping), | |
92 | fYLumping(o.fYLumping), | |
93 | fTotal(0), | |
94 | fEmpty(0), | |
95 | fBasic(0), | |
96 | fEmptyVsTotal(0), | |
97 | fMean(0), | |
98 | fOcc(0), | |
99 | fCorr(0) | |
100 | { | |
101 | Init(); | |
102 | Reset(o.fBasic); | |
103 | } | |
104 | ||
105 | //____________________________________________________________________ | |
106 | AliPoissonCalculator::~AliPoissonCalculator() | |
107 | { | |
108 | // CleanUp(); | |
109 | } | |
110 | ||
111 | //____________________________________________________________________ | |
112 | void | |
113 | AliPoissonCalculator::CleanUp() | |
114 | { | |
115 | if (fTotal) { delete fTotal; fTotal = 0; } | |
116 | if (fEmpty) { delete fEmpty; fEmpty = 0; } | |
117 | if (fBasic) { delete fBasic; fBasic = 0; } | |
118 | if (fEmptyVsTotal) { delete fEmptyVsTotal; fEmptyVsTotal = 0; } | |
119 | if (fMean) { delete fMean; fMean = 0; } | |
120 | if (fOcc) { delete fOcc; fOcc = 0; } | |
121 | if (fCorr) { delete fCorr; fCorr = 0; } | |
122 | } | |
123 | //____________________________________________________________________ | |
124 | AliPoissonCalculator& | |
125 | AliPoissonCalculator::operator=(const AliPoissonCalculator& o) | |
126 | { | |
127 | if (&o == this) return *this; | |
128 | TNamed::operator=(o); | |
129 | fXLumping = o.fXLumping; | |
130 | fYLumping = o.fYLumping; | |
131 | CleanUp(); | |
132 | Init(); | |
133 | Reset(o.fBasic); | |
134 | return *this; | |
135 | } | |
136 | ||
137 | //____________________________________________________________________ | |
138 | void | |
139 | AliPoissonCalculator::Init(Int_t xLumping, Int_t yLumping) | |
140 | { | |
141 | // | |
142 | // Initialize | |
143 | // | |
144 | if (xLumping > 0) fXLumping = xLumping; | |
145 | if (yLumping > 0) fYLumping = yLumping; | |
146 | ||
147 | //Create diagnostics if void | |
148 | if (fEmptyVsTotal) return; | |
149 | ||
150 | MakeOutput(); | |
151 | } | |
152 | //____________________________________________________________________ | |
153 | void | |
154 | AliPoissonCalculator::Define(const TAxis& xaxis, const TAxis& yaxis) | |
155 | { | |
156 | // | |
157 | // Initialize | |
158 | // | |
159 | const Double_t* xBins = xaxis.GetXbins()->GetArray(); | |
160 | const Double_t* yBins = yaxis.GetXbins()->GetArray(); | |
161 | Int_t nX = xaxis.GetNbins(); | |
162 | Int_t nY = yaxis.GetNbins(); | |
163 | Double_t lX = xaxis.GetXmin(); | |
164 | Double_t hX = xaxis.GetXmax(); | |
165 | Double_t lY = yaxis.GetXmin(); | |
166 | Double_t hY = yaxis.GetXmax(); | |
167 | ||
168 | if (xBins) { | |
169 | if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, yBins); | |
170 | else fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, lY, hY); | |
171 | } | |
172 | else { | |
173 | if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, lX, hX, nY, yBins); | |
174 | else fBasic = new TH2D(kBasicN, kBasicT, nX, lX, hX, nY, lY, hY); | |
175 | } | |
176 | fBasic->SetXTitle(xaxis.GetTitle()); | |
177 | fBasic->SetYTitle(yaxis.GetTitle()); | |
178 | ||
179 | Reset(fBasic); | |
180 | } | |
181 | //____________________________________________________________________ | |
182 | void AliPoissonCalculator::MakeOutput() { | |
183 | ||
184 | Int_t n = fXLumping * fYLumping + 1; | |
185 | fEmptyVsTotal = new TH2D("emptyVsTotal", | |
186 | "# of empty # bins vs total # bins", | |
187 | n, -.5, n-.5, n, -.5, n-.5); | |
188 | fEmptyVsTotal->SetXTitle("Total # bins"); | |
189 | fEmptyVsTotal->SetYTitle("# empty bins"); | |
190 | fEmptyVsTotal->SetZTitle("Correlation"); | |
191 | fEmptyVsTotal->SetOption("colz"); | |
192 | fEmptyVsTotal->SetDirectory(0); | |
193 | ||
194 | n = (fXLumping + fYLumping); | |
195 | fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson", | |
196 | 10*n+1, -.1, n+.1); | |
197 | fMean->SetXTitle("#bar{N_{ch}}=log(empty/total)"); | |
198 | fMean->SetYTitle("Events"); | |
199 | fMean->SetFillColor(kRed+1); | |
200 | fMean->SetFillStyle(3001); | |
201 | fMean->SetLineColor(kBlack); | |
202 | fMean->SetDirectory(0); | |
203 | ||
204 | fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)", | |
205 | 101, -.5, 100.5); | |
206 | fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]"); | |
207 | fOcc->SetYTitle("Events"); | |
208 | fOcc->SetFillColor(kBlue+1); | |
209 | fOcc->SetFillStyle(3001); | |
210 | fOcc->SetLineColor(kBlack); | |
211 | fOcc->SetDirectory(0); | |
212 | ||
213 | fCorr = new TH2D("correction", "Correction as function of mean N_{ch}", | |
214 | 10*n+1, -.1, n+.1, 100, 0, 10); | |
215 | fCorr->SetXTitle("#bar{N_{ch}}"); | |
216 | fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})"); | |
217 | fCorr->SetZTitle("Events"); | |
218 | fCorr->SetOption("colz"); | |
219 | fCorr->SetDirectory(0); | |
220 | ||
221 | ||
222 | } | |
223 | //____________________________________________________________________ | |
224 | void | |
225 | AliPoissonCalculator::Output(TList* d) | |
226 | { | |
227 | if (!d) return; | |
228 | if (!fEmptyVsTotal) MakeOutput(); | |
229 | d->Add(fEmptyVsTotal); | |
230 | d->Add(fMean); | |
231 | d->Add(fOcc); | |
232 | d->Add(fCorr); | |
233 | } | |
234 | ||
235 | //____________________________________________________________________ | |
236 | void | |
237 | AliPoissonCalculator::SetLumping(UShort_t nx, UShort_t ny) | |
238 | { | |
239 | if (nx == fXLumping && ny == fYLumping && | |
240 | fEmptyVsTotal && fTotal) | |
241 | // Check if we have something to do. | |
242 | return; | |
243 | CleanUp(); | |
244 | Init(nx, ny); | |
245 | } | |
246 | ||
247 | //____________________________________________________________________ | |
248 | Int_t | |
249 | AliPoissonCalculator::CheckLumping(char which, Int_t nBins, Int_t lumping) const | |
250 | { | |
251 | if ((nBins % lumping) == 0) return lumping; | |
252 | Int_t l = lumping; | |
253 | do { | |
254 | l--; | |
255 | } while (l > 0 && ((nBins % l) != 0)); | |
256 | Warning("CheckLumping", "%c lumping %d is not a divisor of %d, set to %d", | |
257 | which, lumping, nBins, l); | |
258 | return l; | |
259 | } | |
260 | ||
261 | //____________________________________________________________________ | |
262 | void | |
263 | AliPoissonCalculator::Reset(const TH2D* base) | |
264 | { | |
265 | // | |
266 | // Reset histogram | |
267 | // | |
268 | if (fBasic && fTotal && fEmpty) { | |
269 | fBasic->Reset(); | |
270 | fTotal->Reset(); | |
271 | fEmpty->Reset(); | |
272 | return; | |
273 | } | |
274 | ||
275 | if (!base) return; | |
276 | ||
277 | Int_t nXF = base->GetNbinsX(); | |
278 | Double_t xMin = base->GetXaxis()->GetXmin(); | |
279 | Double_t xMax = base->GetXaxis()->GetXmax(); | |
280 | Int_t nYF = base->GetNbinsY(); | |
281 | Double_t yMin = base->GetYaxis()->GetXmin(); | |
282 | Double_t yMax = base->GetYaxis()->GetXmax(); | |
283 | ||
284 | fXLumping = CheckLumping('X', nXF, fXLumping); | |
285 | fYLumping = CheckLumping('Y', nYF, fYLumping); | |
286 | ||
287 | Int_t nY = nYF / fYLumping; | |
288 | Int_t nX = nXF / fXLumping; | |
289 | ||
290 | if (fBasic != base) { | |
291 | fBasic = static_cast<TH2D*>(base->Clone(kBasicN)); | |
292 | fBasic->SetTitle(kBasicT); | |
293 | fBasic->SetDirectory(0); | |
294 | fBasic->Sumw2(); | |
295 | } | |
296 | ||
297 | fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax); | |
298 | fTotal->SetDirectory(0); | |
299 | fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle()); | |
300 | fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle()); | |
301 | fTotal->Sumw2(); | |
302 | ||
303 | fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN)); | |
304 | fEmpty->SetTitle(kEmptyT); | |
305 | fEmpty->SetDirectory(0); | |
306 | // fEmpty->Sumw2(); | |
307 | } | |
308 | ||
309 | //____________________________________________________________________ | |
310 | void | |
311 | AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight) | |
312 | { | |
313 | // | |
314 | // Fill in an observation | |
315 | // | |
316 | // Parameters: | |
317 | // x X value | |
318 | // Y Y value | |
319 | // hit True if hit | |
320 | // weight Weight if this | |
321 | // | |
322 | fTotal->Fill(x, y); | |
323 | if (hit) fBasic->Fill(x, y, weight); | |
324 | else fEmpty->Fill(x, y); | |
325 | } | |
326 | ||
327 | //____________________________________________________________________ | |
328 | Double_t | |
329 | AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const | |
330 | { | |
331 | if (total <= 0) return 0; | |
332 | if (empty < .001) empty = .001; | |
333 | return -TMath::Log(empty/total); | |
334 | } | |
335 | //____________________________________________________________________ | |
336 | Double_t | |
337 | AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const | |
338 | { | |
339 | if (total <= 0) return 0; | |
340 | if (TMath::Abs(empty-total) < .001) empty = total - .001; | |
341 | return 1 / (1 - empty / total); | |
342 | } | |
343 | ||
344 | //____________________________________________________________________ | |
345 | Int_t | |
346 | AliPoissonCalculator::GetReducedXBin(Int_t ix) const | |
347 | { | |
348 | if (!fBasic) return 0; | |
349 | Double_t mx = fBasic->GetXaxis()->GetBinCenter(ix); | |
350 | return GetReducedXBin(mx); | |
351 | } | |
352 | //____________________________________________________________________ | |
353 | Int_t | |
354 | AliPoissonCalculator::GetReducedXBin(Double_t x) const | |
355 | { | |
356 | if (!fEmpty) return 0; | |
357 | return fEmpty->GetXaxis()->FindBin(x); | |
358 | } | |
359 | //____________________________________________________________________ | |
360 | Int_t | |
361 | AliPoissonCalculator::GetReducedYBin(Int_t iy) const | |
362 | { | |
363 | if (!fBasic) return 0; | |
364 | Double_t my = fBasic->GetYaxis()->GetBinCenter(iy); | |
365 | return GetReducedYBin(my); | |
366 | } | |
367 | //____________________________________________________________________ | |
368 | Int_t | |
369 | AliPoissonCalculator::GetReducedYBin(Double_t y) const | |
370 | { | |
371 | if (!fEmpty) return 0; | |
372 | return fEmpty->GetYaxis()->FindBin(y); | |
373 | } | |
374 | ||
375 | ||
376 | ||
377 | //____________________________________________________________________ | |
378 | TH2D* | |
379 | AliPoissonCalculator::Result(Bool_t correct) | |
380 | { | |
381 | // | |
382 | // Calculate result and store in @a output | |
383 | // | |
384 | // Return: | |
385 | // The result histogram (fBase overwritten) | |
386 | // | |
387 | ||
388 | // Double_t total = fXLumping * fYLumping; | |
389 | ||
390 | for (Int_t ix = 1; ix <= fBasic->GetNbinsX(); ix++) { | |
391 | // Double_t x = fBasic->GetXaxis()->GetBinCenter(ix); | |
392 | Int_t jx = GetReducedXBin(ix); // fEmpty->GetXaxis()->FindBin(x); | |
393 | for (Int_t iy = 1; iy <= fBasic->GetNbinsY(); iy++) { | |
394 | // Double_t y = fBasic->GetYaxis()->GetBinCenter(iy); | |
395 | Int_t jy = GetReducedYBin(iy); // fEmpty->GetYaxis()->FindBin(y); | |
396 | Double_t empty = fEmpty->GetBinContent(jx, jy); | |
397 | Double_t total = fTotal->GetBinContent(jx, jy); | |
398 | Double_t hits = fBasic->GetBinContent(ix,iy); | |
399 | // Mean in region of interest | |
400 | Double_t poissonM = CalculateMean(empty, total); | |
401 | Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1); | |
402 | ||
403 | Double_t poissonV = hits * poissonM * poissonC; | |
404 | Double_t poissonE = TMath::Sqrt(poissonV); | |
405 | if(poissonV > 0) poissonE = TMath::Sqrt(poissonV); | |
406 | ||
407 | fBasic->SetBinContent(ix,iy,poissonV); | |
408 | fBasic->SetBinError(ix,iy,poissonE); | |
409 | } | |
410 | } | |
411 | for (Int_t ix = 1; ix <= fEmpty->GetNbinsX(); ix++) { | |
412 | for (Int_t iy = 1; iy <= fEmpty->GetNbinsY(); iy++) { | |
413 | Double_t empty = fEmpty->GetBinContent(ix, iy); | |
414 | Double_t total = fTotal->GetBinContent(ix, iy); | |
415 | Double_t mean = CalculateMean(empty, total); | |
416 | Double_t corr = CalculateCorrection(empty, total); | |
417 | fEmptyVsTotal->Fill(total, empty); | |
418 | fMean->Fill(mean); | |
419 | if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total)); | |
420 | //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean))); | |
421 | fCorr->Fill(mean, corr); | |
422 | } | |
423 | } | |
424 | return fBasic; | |
425 | } | |
426 | ||
427 | //____________________________________________________________________ | |
428 | void | |
429 | AliPoissonCalculator::Print(const Option_t*) const | |
430 | { | |
431 | // | |
432 | // Print information | |
433 | // | |
434 | // Parameters: | |
435 | // option Not used | |
436 | // | |
437 | char ind[gROOT->GetDirLevel()+3]; | |
438 | for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' '; | |
439 | ind[gROOT->GetDirLevel()] = '\0'; | |
440 | std::cout << ind << ClassName() << ": " << GetName() << '\n' | |
441 | << std::boolalpha | |
442 | << ind << " X lumping: " << fXLumping << '\n' | |
443 | << ind << " Y lumping: " << fYLumping << '\n' | |
444 | << std::noboolalpha << std::endl; | |
445 | } | |
446 | //____________________________________________________________________ | |
447 | void | |
448 | AliPoissonCalculator::Browse(TBrowser* b) | |
449 | { | |
450 | // | |
451 | // Browse this object | |
452 | // | |
453 | // Parameters: | |
454 | // b Object to browse | |
455 | // | |
456 | if (fTotal) b->Add(fTotal); | |
457 | if (fEmpty) b->Add(fEmpty); | |
458 | if (fBasic) b->Add(fBasic); | |
459 | if (fEmptyVsTotal) b->Add(fEmptyVsTotal); | |
460 | if (fMean) b->Add(fMean); | |
461 | if (fOcc) b->Add(fOcc); | |
462 | if (fCorr) b->Add(fCorr); | |
463 | } | |
464 | // | |
465 | // EOF | |
466 | // |