AliPhysics  aaf9c62 (aaf9c62)
TestPoisson.C
Go to the documentation of this file.
1 
10 void
12 {
13  if (nBins <= 0) return;
14  if (max <= 0) max = nBins;
15 
16  Double_t dBin = max / nBins;
17  min = - dBin / 2;
18  max = max + dBin / 2;
19  nBins = nBins + 1;
20 }
21 
31 void
32 TestPoisson(Double_t o=.3, bool useWeights=false, bool correct=true)
33 {
34  const char* load = "$ALICE_PHYSICS/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C";
35  if (!gROOT->GetClass("AliAODForwardMult")) {
36  gROOT->Macro(load);
37  gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
38  }
39 
40  // --- Parameters of this script -----------------------------------
41  Int_t nBin = 5; // Our detector matrix size
42  Int_t nMax = TMath::Max(Int_t(nBin * nBin * o + .5)+nBin/2,nBin);
43  Int_t nEv = 10000; // Number of events
44  Double_t mp = o; // The 'hit' probability
45 
46 
47  TH2D* base = new TH2D("base", "Basic histogram",
48  nBin,-.5, nBin-.5, nBin, -.5, nBin-.5);
49  base->SetXTitle("#eta");
50  base->SetYTitle("#varphi");
51  base->SetDirectory(0);
52  base->SetOption("colz");
53 
54  Int_t tN1=nMax; Double_t tMin1; Double_t tMax1;
55  Int_t tN2=nMax*10; Double_t tMin2; Double_t tMax2=nMax;
56  MakeIntegerAxis(tN1, tMin1, tMax1);
57  MakeIntegerAxis(tN2, tMin2, tMax2);
58  TH2D* corr = new TH2D("comp", "Comparison",
59  tN1, tMin1, tMax1, tN2, tMin2, tMax2);
60  corr->SetXTitle("Input");
61  corr->SetYTitle("Poisson");
62  corr->SetDirectory(0);
63  corr->SetOption("colz");
64  corr->SetStats(0);
65  TLine* lcorr = new TLine(0, 0, tMax2, tMax2);
66 
67  Int_t mm = TMath::Max(Int_t(nBin * o + .5),nBin/2);
68  tN2=mm*10; tMax2 = mm;
69  MakeIntegerAxis(tN2, tMin2, tMax2);
70  Info("", "Making mean w/nbins=%d,range=[%f,%f]", tN2, tMin2, tMax2);
71  TH2D* mean = new TH2D("mean", "Mean comparison",
72  tN2, tMin2, tMax2, tN2, tMin2, tMax2);
73  mean->SetXTitle("Input");
74  mean->SetYTitle("Poisson");
75  mean->SetDirectory(0);
76  mean->SetOption("colz");
77  mean->SetStats(0);
78  TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);
79 
80  TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
81  dist->SetXTitle("s");
82  dist->SetYTitle("P(s)");
83  dist->SetFillColor(kRed+1);
84  dist->SetFillStyle(3001);
85  dist->SetDirectory(0);
86 
87  TH1D* diff = new TH1D("diff", "P-T", 100, -25, 25);
88  diff->SetXTitle("Difference");
89  diff->SetFillColor(kRed+1);
90  diff->SetFillStyle(3001);
91  diff->SetYTitle("Prob");
92 
94  c->Init(nBin ,nBin);
95 
96  for (Int_t i = 0; i < nEv; i++) {
97  c->Reset(base);
98  base->Reset();
99 
100  for (Int_t iEta = 0; iEta < nBin; iEta++) {
101  for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
102  // Throw a die
103  Int_t m = gRandom->Poisson(mp);
104  dist->Fill(m);
105 
106  // Fill into our base histogram
107  base->Fill(iEta, iPhi, m);
108 
109  // Fill into poisson calculator
110  c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
111  }
112  }
113  // Calculate the result
114  TH2D* res = c->Result(correct);
115 
116  // Now loop and compare
117  Double_t mBase = 0;
118  Double_t mPois = 0;
119  for (Int_t iEta = 0; iEta < nBin; iEta++) {
120  for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
121  Double_t p = res->GetBinContent(iEta, iPhi);
122  Double_t t = base->GetBinContent(iEta, iPhi);
123 
124  mBase += t;
125  mPois += p;
126  corr->Fill(t, p);
127  diff->Fill(p-t);
128  }
129  }
130  Int_t nn = nBin * nBin;
131  mean->Fill(mBase / nn, mPois / nn);
132  }
133 
134  TCanvas* cc = new TCanvas("c", "c", 900, 900);
135  cc->SetFillColor(0);
136  cc->SetFillStyle(0);
137  cc->SetBorderMode(0);
138  cc->SetRightMargin(0.02);
139  cc->SetTopMargin(0.02);
140  cc->Divide(2,2);
141 
142  TVirtualPad* pp = cc->cd(1);
143  pp->SetFillColor(0);
144  pp->SetFillStyle(0);
145  pp->SetBorderMode(0);
146  pp->SetRightMargin(0.15);
147  pp->SetTopMargin(0.02);
148  pp->SetLogz();
149  pp->SetGridx();
150  pp->SetGridy();
151  corr->Draw();
152  lcorr->Draw();
153 
154  pp = cc->cd(2);
155  pp->SetFillColor(0);
156  pp->SetFillStyle(0);
157  pp->SetBorderMode(0);
158  pp->SetRightMargin(0.02);
159  pp->SetTopMargin(0.02);
160 #if 0
161  c->GetMean()->Draw();
162 #elif 1
163  pp->SetLogy();
164  diff->Draw();
165 #elif 1
166  c->GetOccupancy()->Draw();
167 #else
168  pp->SetLogy();
169  dist->SetStats(0);
170  dist->Scale(1. / dist->Integral());
171  dist->Draw();
172  TH1D* m1 = c->GetMean();
173  m1->Scale(1. / m1->Integral());
174  m1->Draw("same");
175  Double_t eI;
176  Double_t ii = 100 * dist->Integral(2, 0);
177  TLatex* ll = new TLatex(.97, .85,
178  Form("Input #bar{m}: %5.3f", mp));
179  ll->SetNDC();
180  ll->SetTextFont(132);
181  ll->SetTextAlign(31);
182  ll->Draw();
183  ll->DrawLatex(.97, .75, Form("Result #bar{m}: %5.3f", dist->GetMean()));
184  ll->DrawLatex(.97, .65, Form("Occupancy: #int_{1}^{#infty}P(s)ds = %6.2f%%",
185  ii));
186 
187 #endif
188 
189  pp = cc->cd(3);
190  pp->SetFillColor(0);
191  pp->SetFillStyle(0);
192  pp->SetBorderMode(0);
193  pp->SetRightMargin(0.15);
194  pp->SetTopMargin(0.02);
195  pp->SetGridx();
196  pp->SetGridy();
197  c->GetCorrection()->Draw();
198 
199  pp = cc->cd(4);
200  pp->SetFillColor(0);
201  pp->SetFillStyle(0);
202  pp->SetBorderMode(0);
203  pp->SetRightMargin(0.15);
204  pp->SetTopMargin(0.02);
205  pp->SetLogz();
206  pp->SetGridx();
207  pp->SetGridy();
208  mean->Draw();
209  lmean->Draw();
210 
211  cc->cd();
212 }
213 
214 //
215 // EOF
216 //
217 
218 
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
void Reset(const TH2D *base)
TSystem * gSystem
TCanvas * c
Definition: TestFitELoss.C:172
TH2D * Result(Bool_t correct=true)
TH1D * GetOccupancy() const
TRandom * gRandom
void TestPoisson(Double_t o=.3, bool useWeights=false, bool correct=true)
Definition: TestPoisson.C:32
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)
void MakeIntegerAxis(Int_t &nBins, Double_t &min, Double_t &max)
Definition: TestPoisson.C:11
TH2D * GetCorrection() const