PWG2/FORWARD/analysis2/tests/TestPoisson.C

   1 void
   2 MakeIntegerAxis(Int_t& nBins, Double_t& min, Double_t& max)
   3 {
   4   if (nBins <= 0) return;
   5   if (max <= 0)   max = nBins;
   6
   7   Double_t dBin = max / nBins;
   8   min   = - dBin / 2;
   9   max   = max + dBin / 2;
  10   nBins = nBins + 1;
  11 }
  12
  13 void
  14 TestPoisson(Double_t o=.3, bool useWeights=false)
  15 {
  16   const char* load = "$ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/LoadLibs.C";
  17   gROOT->Macro(load);
  18   gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
  19
  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
  25
  26
  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");
  33
  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");
  44   corr->SetStats(0);
  45   TLine* lcorr = new TLine(0, 0, tMax2, tMax2);
  46
  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");
  57   mean->SetStats(0);
  58   TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);
  59
  60   TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
  61   dist->SetXTitle("s");
  62   dist->SetYTitle("P(s)");
  63   dist->SetFillColor(kRed+1);
  64   dist->SetFillStyle(3001);
  65   dist->SetDirectory(0);
  66
  67   AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
  68   c->Init(nBin ,nBin);
  69
  70   for (Int_t i = 0; i < nEv; i++) {
  71     c->Reset(base);
  72     base->Reset();
  73
  74     for (Int_t iEta = 0; iEta < nBin; iEta++) {
  75       for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
  76         // Throw a die
  77         Int_t m = gRandom->Poisson(mp);
  78         dist->Fill(m);
  79
  80         // Fill into our base histogram
  81         base->Fill(iEta, iPhi, m);
  82
  83         // Fill into poisson calculator
  84         c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
  85       }
  86     }
  87     // Calculate the result
  88     TH2D* res = c->Result();
  89
  90     // Now loop and compare
  91     Double_t mBase = 0;
  92     Double_t mPois = 0;
  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);
  97
  98         mBase += t;
  99         mPois += p;
 100         corr->Fill(t, p);
 101       }
 102     }
 103     Int_t nn = nBin * nBin;
 104     mean->Fill(mBase / nn, mPois / nn);
 105   }
 106
 107   TCanvas* cc = new TCanvas("c", "c", 900, 900);
 108   cc->SetFillColor(0);
 109   cc->SetFillStyle(0);
 110   cc->SetBorderMode(0);
 111   cc->SetRightMargin(0.02);
 112   cc->SetTopMargin(0.02);
 113   cc->Divide(2,2);
 114
 115   TVirtualPad* pp = cc->cd(1);
 116   pp->SetFillColor(0);
 117   pp->SetFillStyle(0);
 118   pp->SetBorderMode(0);
 119   pp->SetRightMargin(0.15);
 120   pp->SetTopMargin(0.02);
 121   pp->SetLogz();
 122   pp->SetGridx();
 123   pp->SetGridy();
 124   corr->Draw();
 125   lcorr->Draw();
 126
 127   pp = cc->cd(2);
 128   pp->SetFillColor(0);
 129   pp->SetFillStyle(0);
 130   pp->SetBorderMode(0);
 131   pp->SetRightMargin(0.02);
 132   pp->SetTopMargin(0.02);
 133 #if 1
 134   c->GetMean()->Draw();
 135 #elif 1
 136   c->GetOccupancy()->Draw();
 137 #else
 138   pp->SetLogy();
 139   dist->SetStats(0);
 140   dist->Scale(1. / dist->Integral());
 141   dist->Draw();
 142   TH1D* m1 = c->GetMean();
 143   m1->Scale(1. / m1->Integral());
 144   m1->Draw("same");
 145   Double_t eI;
 146   Double_t ii = 100 * dist->Integral(2, 0);
 147   TLatex* ll = new TLatex(.97, .85,
 148                           Form("Input #bar{m}: %5.3f", mp));
 149   ll->SetNDC();
 150   ll->SetTextFont(132);
 151   ll->SetTextAlign(31);
 152   ll->Draw();
 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%%",
 155                                ii));
 156
 157 #endif
 158
 159   pp = cc->cd(3);
 160   pp->SetFillColor(0);
 161   pp->SetFillStyle(0);
 162   pp->SetBorderMode(0);
 163   pp->SetRightMargin(0.15);
 164   pp->SetTopMargin(0.02);
 165   pp->SetGridx();
 166   pp->SetGridy();
 167   c->GetCorrection()->Draw();
 168
 169   pp = cc->cd(4);
 170   pp->SetFillColor(0);
 171   pp->SetFillStyle(0);
 172   pp->SetBorderMode(0);
 173   pp->SetRightMargin(0.15);
 174   pp->SetTopMargin(0.02);
 175   pp->SetLogz();
 176   pp->SetGridx();
 177   pp->SetGridy();
 178   mean->Draw();
 179   lmean->Draw();
 180
 181   cc->cd();
 182 }
 183
 184 //
 185 // EOF
 186 //
 187
 188