AliPhysics  b69f57c (b69f57c)
AliPoissonCalculator.cxx
Go to the documentation of this file.
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 passed 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 //____________________________________________________________________
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 //____________________________________________________________________
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 //____________________________________________________________________
90  : TNamed(o),
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 //____________________________________________________________________
107 {
108  // CleanUp();
109 }
110
111 //____________________________________________________________________
112 void
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 //____________________________________________________________________
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
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 //____________________________________________________________________
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 void
224 {
225  if (!d) return;
226  if (!fEmptyVsTotal) MakeOutput();
231 }
232
233 //____________________________________________________________________
234 void
236 {
237  if (nx == fXLumping && ny == fYLumping &&
238  fEmptyVsTotal && fTotal)
239  // Check if we have something to do.
240  return;
241  CleanUp();
242  Init(nx, ny);
243 }
244
245 //____________________________________________________________________
246 Int_t
247 AliPoissonCalculator::CheckLumping(char which, Int_t nBins, Int_t lumping) const
248 {
249  if ((nBins % lumping) == 0) return lumping;
250  Int_t l = lumping;
251  do {
252  l--;
253  } while (l > 0 && ((nBins % l) != 0));
254  Warning("CheckLumping", "%c lumping %d is not a divisor of %d, set to %d",
255  which, lumping, nBins, l);
256  return l;
257 }
258
259 //____________________________________________________________________
260 void
262 {
263  //
264  // Reset histogram
265  //
266  if (fBasic && fTotal && fEmpty) {
267  fBasic->Reset();
268  fTotal->Reset();
269  fEmpty->Reset();
270  return;
271  }
272
273  if (!base) return;
274
275  Int_t nXF = base->GetNbinsX();
276  Double_t xMin = base->GetXaxis()->GetXmin();
277  Double_t xMax = base->GetXaxis()->GetXmax();
278  Int_t nYF = base->GetNbinsY();
279  Double_t yMin = base->GetYaxis()->GetXmin();
280  Double_t yMax = base->GetYaxis()->GetXmax();
281
282  fXLumping = CheckLumping('X', nXF, fXLumping);
283  fYLumping = CheckLumping('Y', nYF, fYLumping);
284
285  Int_t nY = nYF / fYLumping;
286  Int_t nX = nXF / fXLumping;
287
288  if (fBasic != base) {
289  fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
290  fBasic->SetTitle(kBasicT);
291  fBasic->SetDirectory(0);
292  fBasic->Sumw2();
293  }
294
295  fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax);
296  fTotal->SetDirectory(0);
297  fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle());
298  fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle());
299  fTotal->Sumw2();
300
301  fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
302  fEmpty->SetTitle(kEmptyT);
303  fEmpty->SetDirectory(0);
304  // fEmpty->Sumw2();
305 }
306
307 //____________________________________________________________________
308 void
310 {
311  //
312  // Fill in an observation
313  //
314  // Parameters:
315  // x X value
316  // Y Y value
317  // hit True if hit
318  // weight Weight if this
319  //
320  fTotal->Fill(x, y);
321  if (hit) fBasic->Fill(x, y, weight);
322  else fEmpty->Fill(x, y);
323 }
324
325 //____________________________________________________________________
326 Double_t
328 {
329  if (total <= 0) return 0;
330  if (empty < .001) empty = .001;
331  return -TMath::Log(empty/total);
332 }
333 //____________________________________________________________________
334 Double_t
336 {
337  if (total <= 0) return 0;
338  if (TMath::Abs(empty-total) < .001) empty = total - .001;
339  return 1 / (1 - empty / total);
340 }
341
342 //____________________________________________________________________
343 Int_t
345 {
346  if (!fBasic) return 0;
347  Double_t mx = fBasic->GetXaxis()->GetBinCenter(ix);
348  return GetReducedXBin(mx);
349 }
350 //____________________________________________________________________
351 Int_t
353 {
354  if (!fEmpty) return 0;
355  return fEmpty->GetXaxis()->FindBin(x);
356 }
357 //____________________________________________________________________
358 Int_t
360 {
361  if (!fBasic) return 0;
362  Double_t my = fBasic->GetYaxis()->GetBinCenter(iy);
363  return GetReducedYBin(my);
364 }
365 //____________________________________________________________________
366 Int_t
368 {
369  if (!fEmpty) return 0;
370  return fEmpty->GetYaxis()->FindBin(y);
371 }
372
373
374
375 //____________________________________________________________________
376 TH2D*
378 {
379  //
380  // Calculate result and store in @a output
381  //
382  // Return:
383  // The result histogram (fBase overwritten)
384  //
385
386  // Double_t total = fXLumping * fYLumping;
387
388  for (Int_t ix = 1; ix <= fBasic->GetNbinsX(); ix++) {
389  // Double_t x = fBasic->GetXaxis()->GetBinCenter(ix);
390  Int_t jx = GetReducedXBin(ix); // fEmpty->GetXaxis()->FindBin(x);
391  for (Int_t iy = 1; iy <= fBasic->GetNbinsY(); iy++) {
392  // Double_t y = fBasic->GetYaxis()->GetBinCenter(iy);
393  Int_t jy = GetReducedYBin(iy); // fEmpty->GetYaxis()->FindBin(y);
394  Double_t empty = fEmpty->GetBinContent(jx, jy);
395  Double_t total = fTotal->GetBinContent(jx, jy);
396  Double_t hits = fBasic->GetBinContent(ix,iy);
397  // Mean in region of interest
398  Double_t poissonM = CalculateMean(empty, total);
399  Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1);
400
401  Double_t poissonV = hits * poissonM * poissonC;
402  Double_t poissonE = TMath::Sqrt(poissonV);
403  if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
404
405  fBasic->SetBinContent(ix,iy,poissonV);
406  fBasic->SetBinError(ix,iy,poissonE);
407  }
408  }
409  return fBasic;
410 }
411
412 //____________________________________________________________________
413 void
415 {
416  for (Int_t ix = 1; ix <= fEmpty->GetNbinsX(); ix++) {
417  for (Int_t iy = 1; iy <= fEmpty->GetNbinsY(); iy++) {
418  Double_t empty = fEmpty->GetBinContent(ix, iy);
419  Double_t total = fTotal->GetBinContent(ix, iy);
420  Double_t mean = CalculateMean(empty, total);
421  Double_t corr = CalculateCorrection(empty, total);
422  fEmptyVsTotal->Fill(total, empty);
423  fMean->Fill(mean);
424  if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total));
425  //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
426  fCorr->Fill(mean, corr);
427  }
428  }
429 }
430 //____________________________________________________________________
431 void
433 {
434  //
435  // Print information
436  //
437  // Parameters:
438  // option Not used
439  //
440  char ind[gROOT->GetDirLevel()+3];
441  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
442  ind[gROOT->GetDirLevel()] = '\0';
443  std::cout << ind << ClassName() << ": " << GetName() << '\n'
444  << std::boolalpha
445  << ind << " X lumping: " << fXLumping << '\n'
446  << ind << " Y lumping: " << fYLumping << '\n'
447  << std::noboolalpha << std::endl;
448 }
449 //____________________________________________________________________
450 void
452 {
453  //
454  // Browse this object
455  //
456  // Parameters:
457  // b Object to browse
458  //
466 }
467 //
468 // EOF
469 //
double Double_t
Definition: External.C:58
Double_t CalculateCorrection(Double_t empty, Double_t total) const
void Reset(const TH2D *base)
void Define(const TAxis &xaxis, const TAxis &yaxis)
void SetLumping(UShort_t nx, UShort_t ny)
TH2D * Result(Bool_t correct=true)
void Print(const Option_t *option="") const
Int_t GetReducedXBin(Int_t ix) const
void Init(Int_t xLumping=-1, Int_t yLumping=-1)
int Int_t
Definition: External.C:63
Definition: External.C:228
Definition: External.C:212
void Fill(UShort_t strip, UShort_t sec, Bool_t hit, Double_t weight=1)
Int_t CheckLumping(char which, Int_t nBins, Int_t lumping) const
Double_t CalculateMean(Double_t empty, Double_t total) const
Int_t GetReducedYBin(Int_t iy) const
unsigned short UShort_t
Definition: External.C:28
AliPoissonCalculator & operator=(const AliPoissonCalculator &o)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Double_t yMin
Double_t yMax