1 #include "AliPoissonCalculator.h"
2 #include "AliForwardCorrectionManager.h"
12 // A class to calculate the multiplicity in @f$(x,y)@f$ bins
13 // using Poisson statistics.
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.
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
24 // \lange m\rangle = -\log\left(\frac{e}{t}\right)
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
30 // c &=& \frac{1}{1 - \exp{-\lange m\rangle}}
31 // &=& \frac{1}{1 - \frac{e}{t}}
33 // and the final number in each cell is then
35 // h_i c \lange m\rangle
37 // where @f$h_i@f$ is the number of hits in the cell @f$i@f$
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";
53 //____________________________________________________________________
54 AliPoissonCalculator::AliPoissonCalculator()
71 //____________________________________________________________________
72 AliPoissonCalculator::AliPoissonCalculator(const char*)
73 : TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
88 //____________________________________________________________________
89 AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
91 fXLumping(o.fXLumping),
92 fYLumping(o.fYLumping),
105 //____________________________________________________________________
106 AliPoissonCalculator::~AliPoissonCalculator()
111 //____________________________________________________________________
113 AliPoissonCalculator::CleanUp()
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; }
123 //____________________________________________________________________
124 AliPoissonCalculator&
125 AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
127 if (&o == this) return *this;
128 TNamed::operator=(o);
129 fXLumping = o.fXLumping;
130 fYLumping = o.fYLumping;
137 //____________________________________________________________________
139 AliPoissonCalculator::Init(Int_t xLumping, Int_t yLumping)
144 if (xLumping > 0) fXLumping = xLumping;
145 if (yLumping > 0) fYLumping = yLumping;
147 //Create diagnostics if void
148 if (fEmptyVsTotal) return;
152 //____________________________________________________________________
154 AliPoissonCalculator::Define(const TAxis& xaxis, const TAxis& yaxis)
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();
169 if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, yBins);
170 else fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, lY, hY);
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);
176 fBasic->SetXTitle(xaxis.GetTitle());
177 fBasic->SetYTitle(yaxis.GetTitle());
181 //____________________________________________________________________
182 void AliPoissonCalculator::MakeOutput() {
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);
194 n = (fXLumping + fYLumping);
195 fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
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);
204 fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
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);
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);
223 //____________________________________________________________________
225 AliPoissonCalculator::Output(TList* d)
228 if (!fEmptyVsTotal) MakeOutput();
229 d->Add(fEmptyVsTotal);
235 //____________________________________________________________________
237 AliPoissonCalculator::SetLumping(UShort_t nx, UShort_t ny)
239 if (nx == fXLumping && ny == fYLumping &&
240 fEmptyVsTotal && fTotal)
241 // Check if we have something to do.
247 //____________________________________________________________________
249 AliPoissonCalculator::CheckLumping(char which, Int_t nBins, Int_t lumping) const
251 if ((nBins % lumping) == 0) return lumping;
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);
261 //____________________________________________________________________
263 AliPoissonCalculator::Reset(const TH2D* base)
268 if (fBasic && fTotal && fEmpty) {
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();
284 fXLumping = CheckLumping('X', nXF, fXLumping);
285 fYLumping = CheckLumping('Y', nYF, fYLumping);
287 Int_t nY = nYF / fYLumping;
288 Int_t nX = nXF / fXLumping;
290 if (fBasic != base) {
291 fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
292 fBasic->SetTitle(kBasicT);
293 fBasic->SetDirectory(0);
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());
303 fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
304 fEmpty->SetTitle(kEmptyT);
305 fEmpty->SetDirectory(0);
309 //____________________________________________________________________
311 AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight)
314 // Fill in an observation
320 // weight Weight if this
323 if (hit) fBasic->Fill(x, y, weight);
324 else fEmpty->Fill(x, y);
327 //____________________________________________________________________
329 AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
331 if (total <= 0) return 0;
332 if (empty < .001) empty = .001;
333 return -TMath::Log(empty/total);
335 //____________________________________________________________________
337 AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
339 if (total <= 0) return 0;
340 if (TMath::Abs(empty-total) < .001) empty = total - .001;
341 return 1 / (1 - empty / total);
344 //____________________________________________________________________
346 AliPoissonCalculator::GetReducedXBin(Int_t ix) const
348 if (!fBasic) return 0;
349 Double_t mx = fBasic->GetXaxis()->GetBinCenter(ix);
350 return GetReducedXBin(mx);
352 //____________________________________________________________________
354 AliPoissonCalculator::GetReducedXBin(Double_t x) const
356 if (!fEmpty) return 0;
357 return fEmpty->GetXaxis()->FindBin(x);
359 //____________________________________________________________________
361 AliPoissonCalculator::GetReducedYBin(Int_t iy) const
363 if (!fBasic) return 0;
364 Double_t my = fBasic->GetYaxis()->GetBinCenter(iy);
365 return GetReducedYBin(my);
367 //____________________________________________________________________
369 AliPoissonCalculator::GetReducedYBin(Double_t y) const
371 if (!fEmpty) return 0;
372 return fEmpty->GetYaxis()->FindBin(y);
377 //____________________________________________________________________
379 AliPoissonCalculator::Result(Bool_t correct)
382 // Calculate result and store in @a output
385 // The result histogram (fBase overwritten)
388 // Double_t total = fXLumping * fYLumping;
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);
403 Double_t poissonV = hits * poissonM * poissonC;
404 Double_t poissonE = TMath::Sqrt(poissonV);
405 if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
407 fBasic->SetBinContent(ix,iy,poissonV);
408 fBasic->SetBinError(ix,iy,poissonE);
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);
419 fOcc->Fill(100 * (1 - empty/total));
420 //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
421 fCorr->Fill(mean, corr);
427 //____________________________________________________________________
429 AliPoissonCalculator::Print(const Option_t*) const
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'
442 << ind << " X lumping: " << fXLumping << '\n'
443 << ind << " Y lumping: " << fYLumping << '\n'
444 << std::noboolalpha << std::endl;
446 //____________________________________________________________________
448 AliPoissonCalculator::Browse(TBrowser* b)
451 // Browse this object
454 // b Object to browse
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);