PWG2/FORWARD/analysis2/tests/TestPoisson.C

   1 /**
   2  *
   3  *
   4  * @param nBins
   5  * @param min
   6  * @param max
   7  *
   8  * @ingroup pwg2_forward_analysis_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 pwg2_forward_analysis_scripts_tests
  29  */
  30 void
  31 TestPoisson(Double_t o=.3, bool useWeights=false)
  32 {
  33   const char* load = "$ALICE_ROOT/PWG2/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   AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
  85   c->Init(nBin ,nBin);
  86
  87   for (Int_t i = 0; i < nEv; i++) {
  88     c->Reset(base);
  89     base->Reset();
  90
  91     for (Int_t iEta = 0; iEta < nBin; iEta++) {
  92       for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
  93         // Throw a die
  94         Int_t m = gRandom->Poisson(mp);
  95         dist->Fill(m);
  96
  97         // Fill into our base histogram
  98         base->Fill(iEta, iPhi, m);
  99
 100         // Fill into poisson calculator
 101         c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
 102       }
 103     }
 104     // Calculate the result
 105     TH2D* res = c->Result();
 106
 107     // Now loop and compare
 108     Double_t mBase = 0;
 109     Double_t mPois = 0;
 110     for (Int_t iEta = 0; iEta < nBin; iEta++) {
 111       for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
 112         Double_t p = res->GetBinContent(iEta, iPhi);
 113         Double_t t = base->GetBinContent(iEta, iPhi);
 114
 115         mBase += t;
 116         mPois += p;
 117         corr->Fill(t, p);
 118       }
 119     }
 120     Int_t nn = nBin * nBin;
 121     mean->Fill(mBase / nn, mPois / nn);
 122   }
 123
 124   TCanvas* cc = new TCanvas("c", "c", 900, 900);
 125   cc->SetFillColor(0);
 126   cc->SetFillStyle(0);
 127   cc->SetBorderMode(0);
 128   cc->SetRightMargin(0.02);
 129   cc->SetTopMargin(0.02);
 130   cc->Divide(2,2);
 131
 132   TVirtualPad* pp = cc->cd(1);
 133   pp->SetFillColor(0);
 134   pp->SetFillStyle(0);
 135   pp->SetBorderMode(0);
 136   pp->SetRightMargin(0.15);
 137   pp->SetTopMargin(0.02);
 138   pp->SetLogz();
 139   pp->SetGridx();
 140   pp->SetGridy();
 141   corr->Draw();
 142   lcorr->Draw();
 143
 144   pp = cc->cd(2);
 145   pp->SetFillColor(0);
 146   pp->SetFillStyle(0);
 147   pp->SetBorderMode(0);
 148   pp->SetRightMargin(0.02);
 149   pp->SetTopMargin(0.02);
 150 #if 1
 151   c->GetMean()->Draw();
 152 #elif 1
 153   c->GetOccupancy()->Draw();
 154 #else
 155   pp->SetLogy();
 156   dist->SetStats(0);
 157   dist->Scale(1. / dist->Integral());
 158   dist->Draw();
 159   TH1D* m1 = c->GetMean();
 160   m1->Scale(1. / m1->Integral());
 161   m1->Draw("same");
 162   Double_t eI;
 163   Double_t ii = 100 * dist->Integral(2, 0);
 164   TLatex* ll = new TLatex(.97, .85,
 165                           Form("Input #bar{m}: %5.3f", mp));
 166   ll->SetNDC();
 167   ll->SetTextFont(132);
 168   ll->SetTextAlign(31);
 169   ll->Draw();
 170   ll->DrawLatex(.97, .75, Form("Result #bar{m}: %5.3f", dist->GetMean()));
 171   ll->DrawLatex(.97, .65, Form("Occupancy: #int_{1}^{#infty}P(s)ds = %6.2f%%",
 172                                ii));
 173
 174 #endif
 175
 176   pp = cc->cd(3);
 177   pp->SetFillColor(0);
 178   pp->SetFillStyle(0);
 179   pp->SetBorderMode(0);
 180   pp->SetRightMargin(0.15);
 181   pp->SetTopMargin(0.02);
 182   pp->SetGridx();
 183   pp->SetGridy();
 184   c->GetCorrection()->Draw();
 185
 186   pp = cc->cd(4);
 187   pp->SetFillColor(0);
 188   pp->SetFillStyle(0);
 189   pp->SetBorderMode(0);
 190   pp->SetRightMargin(0.15);
 191   pp->SetTopMargin(0.02);
 192   pp->SetLogz();
 193   pp->SetGridx();
 194   pp->SetGridy();
 195   mean->Draw();
 196   lmean->Draw();
 197
 198   cc->cd();
 199 }
 200
 201 //
 202 // EOF
 203 //
 204
 205