2 MakeIntegerAxis(Int_t& nBins, Double_t& min, Double_t& max)
4 if (nBins <= 0) return;
5 if (max <= 0) max = nBins;
7 Double_t dBin = max / nBins;
14 TestPoisson(Double_t o=.3, bool useWeights=false)
16 const char* load = "$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C";
18 gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
20 // --- Parameters of this script -----------------------------------
21 Int_t nBin = 5; // Our detector matrix size
22 Int_t nMax = TMath::Max(Int_t(nBin * nBin * o + .5)+nBin/2,nBin);
23 Int_t nEv = 10000; // Number of events
24 Double_t mp = o; // The 'hit' probability
27 TH2D* base = new TH2D("base", "Basic histogram",
28 nBin,-.5, nBin-.5, nBin, -.5, nBin-.5);
29 base->SetXTitle("#eta");
30 base->SetYTitle("#varphi");
31 base->SetDirectory(0);
32 base->SetOption("colz");
34 Int_t tN1=nMax; Double_t tMin1; Double_t tMax1;
35 Int_t tN2=nMax*10; Double_t tMin2; Double_t tMax2=nMax;
36 MakeIntegerAxis(tN1, tMin1, tMax1);
37 MakeIntegerAxis(tN2, tMin2, tMax2);
38 TH2D* corr = new TH2D("comp", "Comparison",
39 tN1, tMin1, tMax1, tN2, tMin2, tMax2);
40 corr->SetXTitle("Input");
41 corr->SetYTitle("Poisson");
42 corr->SetDirectory(0);
43 corr->SetOption("colz");
45 TLine* lcorr = new TLine(0, 0, tMax2, tMax2);
47 Int_t mm = TMath::Max(Int_t(nBin * o + .5),nBin/2);
48 tN2=mm*10; tMax2 = mm;
49 MakeIntegerAxis(tN2, tMin2, tMax2);
50 Info("", "Making mean w/nbins=%d,range=[%f,%f]", tN2, tMin2, tMax2);
51 TH2D* mean = new TH2D("mean", "Mean comparison",
52 tN2, tMin2, tMax2, tN2, tMin2, tMax2);
53 mean->SetXTitle("Input");
54 mean->SetYTitle("Poisson");
55 mean->SetDirectory(0);
56 mean->SetOption("colz");
58 TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);
60 TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
62 dist->SetYTitle("P(s)");
63 dist->SetFillColor(kRed+1);
64 dist->SetFillStyle(3001);
65 dist->SetDirectory(0);
67 AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
70 for (Int_t i = 0; i < nEv; i++) {
74 for (Int_t iEta = 0; iEta < nBin; iEta++) {
75 for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
77 Int_t m = gRandom->Poisson(mp);
80 // Fill into our base histogram
81 base->Fill(iEta, iPhi, m);
83 // Fill into poisson calculator
84 c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
87 // Calculate the result
88 TH2D* res = c->Result();
90 // Now loop and compare
93 for (Int_t iEta = 0; iEta < nBin; iEta++) {
94 for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
95 Double_t p = res->GetBinContent(iEta, iPhi);
96 Double_t t = base->GetBinContent(iEta, iPhi);
103 Int_t nn = nBin * nBin;
104 mean->Fill(mBase / nn, mPois / nn);
107 TCanvas* cc = new TCanvas("c", "c", 900, 900);
110 cc->SetBorderMode(0);
111 cc->SetRightMargin(0.02);
112 cc->SetTopMargin(0.02);
115 TVirtualPad* pp = cc->cd(1);
118 pp->SetBorderMode(0);
119 pp->SetRightMargin(0.15);
120 pp->SetTopMargin(0.02);
130 pp->SetBorderMode(0);
131 pp->SetRightMargin(0.02);
132 pp->SetTopMargin(0.02);
134 c->GetMean()->Draw();
136 c->GetOccupancy()->Draw();
140 dist->Scale(1. / dist->Integral());
142 TH1D* m1 = c->GetMean();
143 m1->Scale(1. / m1->Integral());
146 Double_t ii = 100 * dist->Integral(2, 0);
147 TLatex* ll = new TLatex(.97, .85,
148 Form("Input #bar{m}: %5.3f", mp));
150 ll->SetTextFont(132);
151 ll->SetTextAlign(31);
153 ll->DrawLatex(.97, .75, Form("Result #bar{m}: %5.3f", dist->GetMean()));
154 ll->DrawLatex(.97, .65, Form("Occupancy: #int_{1}^{#infty}P(s)ds = %6.2f%%",
162 pp->SetBorderMode(0);
163 pp->SetRightMargin(0.15);
164 pp->SetTopMargin(0.02);
167 c->GetCorrection()->Draw();
172 pp->SetBorderMode(0);
173 pp->SetRightMargin(0.15);
174 pp->SetTopMargin(0.02);