PWGLF/FORWARD/analysis2/tests/TestPoisson.C

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