PWG2/FORWARD/analysis2/AliPoissonCalculator.cxx

   1 #include "AliPoissonCalculator.h"
   2 #include <TH2D.h>
   3 #include <TBrowser.h>
   4 #include <TROOT.h>
   5 #include <TMath.h>
   6 #include <TList.h>
   7 #include <iostream>
   8
   9 //
  10 // A class to calculate the multiplicity in @f$(\eta,\varphi)@f$ bins
  11 // using Poisson statistics.
  12 //
  13 // The input is assumed to be binned in @f$(\eta,\varphi)@f$ as
  14 // described by the 2D histogram passwd to the Reset member function.
  15 //
  16 // The data is grouped in to regions as defined by the parameters
  17 // fEtaLumping and fPhiLumping.  The total number of cells and number
  18 // of empty cells is then calculate in each region.  The mean
  19 // multiplicity over the region is then determined as
  20 //
  21 // @f[
  22 // \lange m\rangle = -\log\left(\frac{e}{t}\right)
  23 // @f]
  24 // where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
  25 // total number of cells in the region.  A correction for counting
  26 // statistics, is then applied
  27 // @f{eqnarray*}
  28 //    c &=& \frac{1}{1 - \exp{-\lange m\rangle}}
  29 //      &=& \frac{1}{1 - \frac{e}{t}}
  30 // @f{eqnarray*}
  31 // and the final number in each cell is then
  32 // @f[
  33 //   h_i c \lange m\rangle
  34 // @f]
  35 // where @f$h_i@f$ is the number of hits in the cell @f$i@f$
  36 //
  37 //
  38
  39 //____________________________________________________________________
  40 AliPoissonCalculator::AliPoissonCalculator()
  41   : TNamed(),
  42     fEtaLumping(5),
  43     fPhiLumping(5),
  44     fTotal(0),
  45     fEmpty(0),
  46     fBasic(0),
  47     fEmptyVsTotal(0),
  48     fMean(0),
  49     fOcc(0),
  50     fCorr(0)
  51 {
  52   //
  53   // CTOR
  54   //
  55 }
  56
  57 //____________________________________________________________________
  58 AliPoissonCalculator::AliPoissonCalculator(const char*)
  59   : TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
  60     fEtaLumping(5),
  61     fPhiLumping(5),
  62     fTotal(0),
  63     fEmpty(0),
  64     fBasic(0),
  65     fEmptyVsTotal(0),
  66     fMean(0),
  67     fOcc(0),
  68     fCorr(0)
  69 {
  70   //
  71   // CTOR
  72   //
  73 }
  74 //____________________________________________________________________
  75 AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
  76   : TNamed(o),
  77     fEtaLumping(o.fEtaLumping),
  78     fPhiLumping(o.fPhiLumping),
  79     fTotal(0),
  80     fEmpty(0),
  81     fBasic(0),
  82     fEmptyVsTotal(0),
  83     fMean(0),
  84     fOcc(0),
  85     fCorr(0)
  86 {
  87   Init();
  88   Reset(o.fBasic);
  89 }
  90
  91 //____________________________________________________________________
  92 AliPoissonCalculator::~AliPoissonCalculator()
  93 {
  94   CleanUp();
  95 }
  96
  97 //____________________________________________________________________
  98 void
  99 AliPoissonCalculator::CleanUp()
 100 {
 101   if (fTotal)        delete fTotal;        fTotal        = 0;
 102   if (fEmpty)        delete fEmpty;        fEmpty        = 0;
 103   if (fBasic)        delete fBasic;        fBasic        = 0;
 104   if (fEmptyVsTotal) delete fEmptyVsTotal; fEmptyVsTotal = 0;
 105   if (fMean)         delete fMean;         fMean         = 0;
 106   if (fOcc)          delete fOcc;          fOcc          = 0;
 107   if (fCorr)         delete fCorr;         fCorr         = 0;
 108 }
 109 //____________________________________________________________________
 110 AliPoissonCalculator&
 111 AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
 112 {
 113   TNamed::operator=(o);
 114   fEtaLumping = o.fEtaLumping;
 115   fPhiLumping = o.fPhiLumping;
 116   CleanUp();
 117   Init(-1,-1);
 118   Reset(o.fBasic);
 119   return *this;
 120 }
 121
 122 //____________________________________________________________________
 123 void
 124 AliPoissonCalculator::Init(Int_t etaLumping, Int_t phiLumping)
 125 {
 126   //
 127   // Initialize
 128   //
 129   if (etaLumping > 0) SetEtaLumping(etaLumping);
 130   if (phiLumping > 0) SetPhiLumping(phiLumping);
 131   if (fEmptyVsTotal) return;
 132
 133   Int_t n = fEtaLumping * fPhiLumping + 1;
 134   fEmptyVsTotal = new TH2D("emptyVsTotal",
 135                            "# of empty # bins vs total # bins",
 136                            n, -.5, n-.5, n, -.5, n-.5);
 137   fEmptyVsTotal->SetXTitle("Total # bins");
 138   fEmptyVsTotal->SetYTitle("# empty bins");
 139   fEmptyVsTotal->SetZTitle("Correlation");
 140   fEmptyVsTotal->SetOption("colz");
 141   fEmptyVsTotal->SetDirectory(0);
 142
 143   n = (fEtaLumping + fPhiLumping);
 144   fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
 145                    10*n+1, -.1, n+.1);
 146   fMean->SetXTitle("#bar{N_{ch}}=#log(empty/total)");
 147   fMean->SetYTitle("Events");
 148   fMean->SetFillColor(kRed+1);
 149   fMean->SetFillStyle(3001);
 150   fMean->SetLineColor(kBlack);
 151   fMean->SetDirectory(0);
 152
 153   fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
 154                   100, 0, 100);
 155   fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
 156   fOcc->SetYTitle("Events");
 157   fOcc->SetFillColor(kBlue+1);
 158   fOcc->SetFillStyle(3001);
 159   fOcc->SetLineColor(kBlack);
 160   fOcc->SetDirectory(0);
 161
 162   fCorr = new TH2D("correction", "Correction as function of mean N_{ch}",
 163                    10*n+1, -.1, n+.1, 100, 0, 10);
 164   fCorr->SetXTitle("#bar{N_{ch}}");
 165   fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})");
 166   fCorr->SetZTitle("Events");
 167   fCorr->SetOption("colz");
 168   fCorr->SetDirectory(0);
 169 }
 170
 171 //____________________________________________________________________
 172 void
 173 AliPoissonCalculator::Output(TList* d)
 174 {
 175   if (!d) return;
 176   if (!fEmptyVsTotal) Init();
 177   d->Add(fEmptyVsTotal);
 178   d->Add(fMean);
 179   d->Add(fOcc);
 180   d->Add(fCorr);
 181 }
 182
 183 //____________________________________________________________________
 184 void
 185 AliPoissonCalculator::Reset(const TH2D* base)
 186 {
 187   //
 188   // Reset histogram
 189   //
 190   if (!base) return;
 191   if (fBasic && fTotal && fEmpty) {
 192     fBasic->Reset();
 193     fTotal->Reset();
 194     fEmpty->Reset();
 195     return;
 196   }
 197
 198   Int_t    nEtaF   = base->GetNbinsX();
 199   Int_t    nEta    = nEtaF / fEtaLumping;
 200   Double_t etaMin  = base->GetXaxis()->GetXmin();
 201   Double_t etaMax  = base->GetXaxis()->GetXmax();
 202   Int_t    nPhiF   = base->GetNbinsY();
 203   Int_t    nPhi    = nPhiF / fPhiLumping;
 204   Double_t phiMin  = base->GetYaxis()->GetXmin();
 205   Double_t phiMax  = base->GetYaxis()->GetXmax();
 206
 207   fTotal = new TH2D("total", "Total number of bins/region",
 208                     nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
 209   fEmpty = new TH2D("empty", "Empty number of bins/region",
 210                     nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
 211   fBasic = new TH2D("basic", "Basic number of hits",
 212                     nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
 213
 214   fTotal->SetDirectory(0);
 215   fEmpty->SetDirectory(0);
 216   fBasic->SetDirectory(0);
 217   fTotal->SetXTitle("#eta");
 218   fEmpty->SetXTitle("#eta");
 219   fBasic->SetXTitle("#eta");
 220   fTotal->SetYTitle("#varphi [radians]");
 221   fEmpty->SetYTitle("#varphi [radians]");
 222   fBasic->SetYTitle("#varphi [radians]");
 223   fTotal->Sumw2();
 224   fEmpty->Sumw2();
 225   fBasic->Sumw2();
 226 }
 227
 228 //____________________________________________________________________
 229 void
 230 AliPoissonCalculator::Fill(Double_t eta, Double_t phi, Bool_t hit,
 231                            Double_t weight)
 232 {
 233   //
 234   // Fill in an observation
 235   //
 236   // Parameters:
 237   //    eta     Eta value
 238   //    phi     Phi value
 239   //    hit     True if hit
 240   //    weight  Weight if this
 241   //
 242   fTotal->Fill(eta, phi);
 243   if (hit) fBasic->Fill(eta, phi, weight);
 244   else     fEmpty->Fill(eta, phi);
 245 }
 246
 247 //____________________________________________________________________
 248 Double_t
 249 AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
 250 {
 251   if (total <= 0) return 0;
 252   if (empty < .001) empty = .001;
 253   return -TMath::Log(empty/total);
 254 }
 255 //____________________________________________________________________
 256 Double_t
 257 AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
 258 {
 259   if (total <= 0) return 0;
 260   if (TMath::Abs(empty-total) < .001) empty = total - .001;
 261   return 1 / (1 - empty / total);
 262 }
 263
 264 //____________________________________________________________________
 265 TH2D*
 266 AliPoissonCalculator::Result()
 267 {
 268   //
 269   // Calculate result and store in @a output
 270   //
 271   // Return:
 272   //    The result histogram (fBase overwritten)
 273   //
 274   for (Int_t ieta = 1; ieta <= fBasic->GetNbinsX(); ieta++) {
 275     for (Int_t iphi = 1; iphi <= fBasic->GetNbinsY(); iphi++) {
 276       Double_t eta      = fBasic->GetXaxis()->GetBinCenter(ieta);
 277       Double_t phi      = fBasic->GetYaxis()->GetBinCenter(iphi);
 278       Int_t    jeta     = fEmpty->GetXaxis()->FindBin(eta);
 279       Int_t    jphi     = fEmpty->GetYaxis()->FindBin(phi);
 280       Double_t empty    = fEmpty->GetBinContent(jeta, jphi);
 281       Double_t total    = fTotal->GetBinContent(jeta, jphi);
 282       Double_t hits     = fBasic->GetBinContent(ieta,iphi);
 283
 284       // Mean in region of interest
 285       Double_t poissonM = CalculateMean(empty, total);
 286       Double_t poissonC = CalculateCorrection(empty, total);
 287       Double_t poissonV = hits * poissonM * poissonC;
 288       Double_t poissonE = TMath::Sqrt(poissonV);
 289       if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
 290
 291       fBasic->SetBinContent(ieta,iphi,poissonV);
 292       fBasic->SetBinError(ieta,iphi,poissonE);
 293     }
 294   }
 295   for (Int_t ieta = 1; ieta <= fEmpty->GetNbinsX(); ieta++) {
 296     for (Int_t iphi = 1; iphi <= fEmpty->GetNbinsY(); iphi++) {
 297       Double_t empty    = fEmpty->GetBinContent(ieta, iphi);
 298       Double_t total    = fTotal->GetBinContent(ieta, iphi);
 299       Double_t mean     = CalculateMean(empty, total);
 300       Double_t corr     = CalculateCorrection(empty, total);
 301       fEmptyVsTotal->Fill(total, empty);
 302       fMean->Fill(mean);
 303       fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
 304       fCorr->Fill(mean, corr);
 305     }
 306   }
 307   return fBasic;
 308 }
 309
 310 //____________________________________________________________________
 311 void
 312 AliPoissonCalculator::Print(const Option_t*) const
 313 {
 314   //
 315   // Print information
 316   //
 317   // Parameters:
 318   //    option Not used
 319   //
 320   char ind[gROOT->GetDirLevel()+3];
 321   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 322   ind[gROOT->GetDirLevel()] = '\0';
 323   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 324             << std::boolalpha
 325             << ind << " Eta lumping:            " << fEtaLumping << '\n'
 326             << ind << " Phi lumping:            " << fPhiLumping << '\n'
 327             << std::noboolalpha << std::endl;
 328 }
 329 //____________________________________________________________________
 330 void
 331 AliPoissonCalculator::Browse(TBrowser* b)
 332 {
 333   //
 334   // Browse this object
 335   //
 336   // Parameters:
 337   //    b Object to browse
 338   //
 339   if (fTotal)        b->Add(fTotal);
 340   if (fEmpty)        b->Add(fEmpty);
 341   if (fBasic)        b->Add(fBasic);
 342   if (fEmptyVsTotal) b->Add(fEmptyVsTotal);
 343   if (fMean)         b->Add(fMean);
 344   if (fOcc)          b->Add(fOcc);
 345   if (fCorr)         b->Add(fCorr);
 346 }
 347 //
 348 // EOF
 349 //