PWGLF/FORWARD/analysis2/AliPoissonCalculator.cxx

   1 #include "AliPoissonCalculator.h"
   2 // #include "AliForwardCorrectionManager.h"
   3 #include <TH2D.h>
   4 #include <TBrowser.h>
   5 #include <TROOT.h>
   6 #include <TMath.h>
   7 #include <TList.h>
   8 #include <iostream>
   9 #include <TAxis.h>
  10
  11 //
  12 // A class to calculate the multiplicity in @f$(x,y)@f$ bins
  13 // using Poisson statistics.
  14 //
  15 // The input is assumed to be binned in @f$(x,y)@f$ as described by
  16 // the 2D histogram passed to the Reset member function.
  17 //
  18 // The data is grouped in to regions as defined by the parameters
  19 // fXLumping and fYLumping.  The total number of cells and number of
  20 // empty cells is then calculate in each region.  The mean
  21 // multiplicity over the region is then determined as
  22 //
  23 // @f[
  24 // \lange m\rangle = -\log\left(\frac{e}{t}\right)
  25 // @f]
  26 // where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
  27 // total number of cells in the region.  A correction for counting
  28 // statistics, is then applied
  29 // @f{eqnarray*}
  30 //    c &=& \frac{1}{1 - \exp{-\lange m\rangle}}
  31 //      &=& \frac{1}{1 - \frac{e}{t}}
  32 // @f{eqnarray*}
  33 // and the final number in each cell is then
  34 // @f[
  35 //   h_i c \lange m\rangle
  36 // @f]
  37 // where @f$h_i@f$ is the number of hits in the cell @f$i@f$
  38 //
  39 //
  40
  41 namespace {
  42   const char* kBasicN = "basic";
  43   const char* kEmptyN = "empty";
  44   const char* kTotalN = "total";
  45   const char* kBasicT = "Basic number of hits";
  46   const char* kEmptyT = "Empty number of bins/region";
  47   const char* kTotalT = "Total number of bins/region";
  48 }
  49
  50
  51
  52
  53 //____________________________________________________________________
  54 AliPoissonCalculator::AliPoissonCalculator()
  55   : TNamed(),
  56     fXLumping(32),
  57     fYLumping(4),
  58     fTotal(0),
  59     fEmpty(0),
  60     fBasic(0),
  61     fEmptyVsTotal(0),
  62     fMean(0),
  63     fOcc(0),
  64     fCorr(0)
  65 {
  66   //
  67   // CTOR
  68   //
  69 }
  70
  71 //____________________________________________________________________
  72 AliPoissonCalculator::AliPoissonCalculator(const char*)
  73   : TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
  74     fXLumping(32),
  75     fYLumping(4),
  76     fTotal(0),
  77     fEmpty(0),
  78     fBasic(0),
  79     fEmptyVsTotal(0),
  80     fMean(0),
  81     fOcc(0),
  82     fCorr(0)
  83 {
  84   //
  85   // CTOR
  86   //
  87 }
  88 //____________________________________________________________________
  89 AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
  90   : TNamed(o),
  91     fXLumping(o.fXLumping),
  92     fYLumping(o.fYLumping),
  93     fTotal(0),
  94     fEmpty(0),
  95     fBasic(0),
  96     fEmptyVsTotal(0),
  97     fMean(0),
  98     fOcc(0),
  99     fCorr(0)
 100 {
 101   Init();
 102   Reset(o.fBasic);
 103 }
 104
 105 //____________________________________________________________________
 106 AliPoissonCalculator::~AliPoissonCalculator()
 107 {
 108   // CleanUp();
 109 }
 110
 111 //____________________________________________________________________
 112 void
 113 AliPoissonCalculator::CleanUp()
 114 {
 115   if (fTotal)        { delete fTotal;        fTotal        = 0; }
 116   if (fEmpty)        { delete fEmpty;        fEmpty        = 0; }
 117   if (fBasic)        { delete fBasic;        fBasic        = 0; }
 118   if (fEmptyVsTotal) { delete fEmptyVsTotal; fEmptyVsTotal = 0; }
 119   if (fMean)         { delete fMean;         fMean         = 0; }
 120   if (fOcc)          { delete fOcc;          fOcc          = 0; }
 121   if (fCorr)         { delete fCorr;         fCorr         = 0; }
 122 }
 123 //____________________________________________________________________
 124 AliPoissonCalculator&
 125 AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
 126 {
 127   if (&o == this) return *this;
 128   TNamed::operator=(o);
 129   fXLumping = o.fXLumping;
 130   fYLumping = o.fYLumping;
 131   CleanUp();
 132   Init();
 133   Reset(o.fBasic);
 134   return *this;
 135 }
 136
 137 //____________________________________________________________________
 138 void
 139 AliPoissonCalculator::Init(Int_t xLumping, Int_t yLumping)
 140 {
 141   //
 142   // Initialize
 143   //
 144   if (xLumping > 0) fXLumping = xLumping;
 145   if (yLumping > 0) fYLumping = yLumping;
 146
 147   //Create diagnostics if void
 148   if (fEmptyVsTotal) return;
 149
 150   MakeOutput();
 151 }
 152 //____________________________________________________________________
 153 void
 154 AliPoissonCalculator::Define(const TAxis& xaxis, const TAxis& yaxis)
 155 {
 156   //
 157   // Initialize
 158   //
 159   const Double_t* xBins = xaxis.GetXbins()->GetArray();
 160   const Double_t* yBins = yaxis.GetXbins()->GetArray();
 161   Int_t           nX    = xaxis.GetNbins();
 162   Int_t           nY    = yaxis.GetNbins();
 163   Double_t        lX    = xaxis.GetXmin();
 164   Double_t        hX    = xaxis.GetXmax();
 165   Double_t        lY    = yaxis.GetXmin();
 166   Double_t        hY    = yaxis.GetXmax();
 167
 168   if (xBins) {
 169     if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, yBins);
 170     else       fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, lY, hY);
 171   }
 172   else {
 173     if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, lX, hX, nY, yBins);
 174     else       fBasic = new TH2D(kBasicN, kBasicT, nX, lX, hX, nY, lY, hY);
 175   }
 176   fBasic->SetXTitle(xaxis.GetTitle());
 177   fBasic->SetYTitle(yaxis.GetTitle());
 178
 179   Reset(fBasic);
 180 }
 181 //____________________________________________________________________
 182 void AliPoissonCalculator::MakeOutput()
 183 {
 184   Int_t n = fXLumping * fYLumping + 1;
 185   fEmptyVsTotal = new TH2D("emptyVsTotal",
 186                            "# of empty # bins vs total # bins",
 187                            n, -.5, n-.5, n, -.5, n-.5);
 188   fEmptyVsTotal->SetXTitle("Total # bins");
 189   fEmptyVsTotal->SetYTitle("# empty bins");
 190   fEmptyVsTotal->SetZTitle("Correlation");
 191   fEmptyVsTotal->SetOption("colz");
 192   fEmptyVsTotal->SetDirectory(0);
 193
 194   n = (fXLumping + fYLumping);
 195   fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
 196                    10*n+1, -.1, n+.1);
 197   fMean->SetXTitle("#bar{N_{ch}}=log(empty/total)");
 198   fMean->SetYTitle("Events");
 199   fMean->SetFillColor(kRed+1);
 200   fMean->SetFillStyle(3001);
 201   fMean->SetLineColor(kBlack);
 202   fMean->SetDirectory(0);
 203
 204   fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
 205                   101, -.5, 100.5);
 206   fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
 207   fOcc->SetYTitle("Events");
 208   fOcc->SetFillColor(kBlue+1);
 209   fOcc->SetFillStyle(3001);
 210   fOcc->SetLineColor(kBlack);
 211   fOcc->SetDirectory(0);
 212
 213   fCorr = new TH2D("correction", "Correction as function of mean N_{ch}",
 214                    10*n+1, -.1, n+.1, 100, 0, 10);
 215   fCorr->SetXTitle("#bar{N_{ch}}");
 216   fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})");
 217   fCorr->SetZTitle("Events");
 218   fCorr->SetOption("colz");
 219   fCorr->SetDirectory(0);
 220 }
 221 //____________________________________________________________________
 222 void
 223 AliPoissonCalculator::Output(TList* d)
 224 {
 225   if (!d) return;
 226   if (!fEmptyVsTotal) MakeOutput();
 227   d->Add(fEmptyVsTotal);
 228   d->Add(fMean);
 229   d->Add(fOcc);
 230   d->Add(fCorr);
 231 }
 232
 233 //____________________________________________________________________
 234 void
 235 AliPoissonCalculator::SetLumping(UShort_t nx, UShort_t ny)
 236 {
 237   if (nx == fXLumping && ny == fYLumping &&
 238       fEmptyVsTotal && fTotal)
 239     // Check if we have something to do.
 240     return;
 241   CleanUp();
 242   Init(nx, ny);
 243 }
 244
 245 //____________________________________________________________________
 246 Int_t
 247 AliPoissonCalculator::CheckLumping(char which, Int_t nBins, Int_t lumping) const
 248 {
 249   if ((nBins % lumping) == 0) return lumping;
 250   Int_t l = lumping;
 251   do {
 252     l--;
 253   } while (l > 0 && ((nBins % l) != 0));
 254   Warning("CheckLumping", "%c lumping %d is not a divisor of %d, set to %d",
 255           which, lumping, nBins, l);
 256   return l;
 257 }
 258
 259 //____________________________________________________________________
 260 void
 261 AliPoissonCalculator::Reset(const TH2D* base)
 262 {
 263   //
 264   // Reset histogram
 265   //
 266   if (fBasic && fTotal && fEmpty) {
 267     fBasic->Reset();
 268     fTotal->Reset();
 269     fEmpty->Reset();
 270     return;
 271   }
 272
 273   if (!base) return;
 274
 275   Int_t    nXF   = base->GetNbinsX();
 276   Double_t xMin  = base->GetXaxis()->GetXmin();
 277   Double_t xMax  = base->GetXaxis()->GetXmax();
 278   Int_t    nYF   = base->GetNbinsY();
 279   Double_t yMin  = base->GetYaxis()->GetXmin();
 280   Double_t yMax  = base->GetYaxis()->GetXmax();
 281
 282   fXLumping = CheckLumping('X', nXF, fXLumping);
 283   fYLumping = CheckLumping('Y', nYF, fYLumping);
 284
 285   Int_t    nY    = nYF / fYLumping;
 286   Int_t    nX    = nXF / fXLumping;
 287
 288   if (fBasic != base) {
 289     fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
 290     fBasic->SetTitle(kBasicT);
 291     fBasic->SetDirectory(0);
 292     fBasic->Sumw2();
 293   }
 294
 295   fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax);
 296   fTotal->SetDirectory(0);
 297   fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle());
 298   fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle());
 299   fTotal->Sumw2();
 300
 301   fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
 302   fEmpty->SetTitle(kEmptyT);
 303   fEmpty->SetDirectory(0);
 304   // fEmpty->Sumw2();
 305 }
 306
 307 //____________________________________________________________________
 308 void
 309 AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight)
 310 {
 311   //
 312   // Fill in an observation
 313   //
 314   // Parameters:
 315   //    x       X value
 316   //    Y       Y value
 317   //    hit     True if hit
 318   //    weight  Weight if this
 319   //
 320   fTotal->Fill(x, y);
 321   if (hit) fBasic->Fill(x, y, weight);
 322   else     fEmpty->Fill(x, y);
 323 }
 324
 325 //____________________________________________________________________
 326 Double_t
 327 AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
 328 {
 329   if (total <= 0) return 0;
 330   if (empty < .001) empty = .001;
 331   return -TMath::Log(empty/total);
 332 }
 333 //____________________________________________________________________
 334 Double_t
 335 AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
 336 {
 337   if (total <= 0) return 0;
 338   if (TMath::Abs(empty-total) < .001) empty = total - .001;
 339   return 1 / (1 - empty / total);
 340 }
 341
 342 //____________________________________________________________________
 343 Int_t
 344 AliPoissonCalculator::GetReducedXBin(Int_t ix) const
 345 {
 346   if (!fBasic) return 0;
 347   Double_t mx = fBasic->GetXaxis()->GetBinCenter(ix);
 348   return GetReducedXBin(mx);
 349 }
 350 //____________________________________________________________________
 351 Int_t
 352 AliPoissonCalculator::GetReducedXBin(Double_t x) const
 353 {
 354   if (!fEmpty) return 0;
 355   return fEmpty->GetXaxis()->FindBin(x);
 356 }
 357 //____________________________________________________________________
 358 Int_t
 359 AliPoissonCalculator::GetReducedYBin(Int_t iy) const
 360 {
 361   if (!fBasic) return 0;
 362   Double_t my = fBasic->GetYaxis()->GetBinCenter(iy);
 363   return GetReducedYBin(my);
 364 }
 365 //____________________________________________________________________
 366 Int_t
 367 AliPoissonCalculator::GetReducedYBin(Double_t y) const
 368 {
 369   if (!fEmpty) return 0;
 370   return fEmpty->GetYaxis()->FindBin(y);
 371 }
 372
 373
 374
 375 //____________________________________________________________________
 376 TH2D*
 377 AliPoissonCalculator::Result(Bool_t correct)
 378 {
 379   //
 380   // Calculate result and store in @a output
 381   //
 382   // Return:
 383   //    The result histogram (fBase overwritten)
 384   //
 385
 386   // Double_t total = fXLumping * fYLumping;
 387
 388   for (Int_t ix = 1; ix <= fBasic->GetNbinsX(); ix++) {
 389     // Double_t x        = fBasic->GetXaxis()->GetBinCenter(ix);
 390     Int_t    jx       = GetReducedXBin(ix); // fEmpty->GetXaxis()->FindBin(x);
 391     for (Int_t iy = 1; iy <= fBasic->GetNbinsY(); iy++) {
 392       // Double_t y        = fBasic->GetYaxis()->GetBinCenter(iy);
 393       Int_t    jy       = GetReducedYBin(iy); // fEmpty->GetYaxis()->FindBin(y);
 394       Double_t empty    = fEmpty->GetBinContent(jx, jy);
 395       Double_t total    = fTotal->GetBinContent(jx, jy);
 396       Double_t hits     = fBasic->GetBinContent(ix,iy);
 397       // Mean in region of interest
 398       Double_t poissonM = CalculateMean(empty, total);
 399       Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1);
 400
 401       Double_t poissonV = hits * poissonM * poissonC;
 402       Double_t poissonE = TMath::Sqrt(poissonV);
 403       if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
 404
 405       fBasic->SetBinContent(ix,iy,poissonV);
 406       fBasic->SetBinError(ix,iy,poissonE);
 407     }
 408   }
 409   return fBasic;
 410 }
 411
 412 //____________________________________________________________________
 413 void
 414 AliPoissonCalculator::FillDiagnostics()
 415 {
 416   for (Int_t ix = 1; ix <= fEmpty->GetNbinsX(); ix++) {
 417     for (Int_t iy = 1; iy <= fEmpty->GetNbinsY(); iy++) {
 418       Double_t empty    = fEmpty->GetBinContent(ix, iy);
 419       Double_t total    = fTotal->GetBinContent(ix, iy);
 420       Double_t mean     = CalculateMean(empty, total);
 421       Double_t corr     = CalculateCorrection(empty, total);
 422       fEmptyVsTotal->Fill(total, empty);
 423       fMean->Fill(mean);
 424       if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total));
 425       //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
 426       fCorr->Fill(mean, corr);
 427     }
 428   }
 429 }
 430 //____________________________________________________________________
 431 void
 432 AliPoissonCalculator::Print(const Option_t*) const
 433 {
 434   //
 435   // Print information
 436   //
 437   // Parameters:
 438   //    option Not used
 439   //
 440   char ind[gROOT->GetDirLevel()+3];
 441   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 442   ind[gROOT->GetDirLevel()] = '\0';
 443   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 444             << std::boolalpha
 445             << ind << " X lumping:              " << fXLumping << '\n'
 446             << ind << " Y lumping:              " << fYLumping << '\n'
 447             << std::noboolalpha << std::endl;
 448 }
 449 //____________________________________________________________________
 450 void
 451 AliPoissonCalculator::Browse(TBrowser* b)
 452 {
 453   //
 454   // Browse this object
 455   //
 456   // Parameters:
 457   //    b Object to browse
 458   //
 459   if (fTotal)        b->Add(fTotal);
 460   if (fEmpty)        b->Add(fEmpty);
 461   if (fBasic)        b->Add(fBasic);
 462   if (fEmptyVsTotal) b->Add(fEmptyVsTotal);
 463   if (fMean)         b->Add(fMean);
 464   if (fOcc)          b->Add(fOcc);
 465   if (fCorr)         b->Add(fCorr);
 466 }
 467 //
 468 // EOF
 469 //