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 passwd 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                   1000, 0, 100);
 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 }
 223 //____________________________________________________________________
 224 void
 225 AliPoissonCalculator::Output(TList* d)
 226 {
 227   if (!d) return;
 228   if (!fEmptyVsTotal) MakeOutput();
 229   d->Add(fEmptyVsTotal);
 230   d->Add(fMean);
 231   d->Add(fOcc);
 232   d->Add(fCorr);
 233 }
 234
 235 //____________________________________________________________________
 236 void
 237 AliPoissonCalculator::SetLumping(UShort_t nx, UShort_t ny)
 238 {
 239   if (nx == fXLumping && ny == fYLumping &&
 240       fEmptyVsTotal && fTotal)
 241     // Check if we have something to do.
 242     return;
 243   CleanUp();
 244   Init(nx, ny);
 245 }
 246
 247 //____________________________________________________________________
 248 Int_t
 249 AliPoissonCalculator::CheckLumping(char which, Int_t nBins, Int_t lumping) const
 250 {
 251   if ((nBins % lumping) == 0) return lumping;
 252   Int_t l = lumping;
 253   do {
 254     l--;
 255   } while (l > 0 && ((nBins % l) != 0));
 256   Warning("CheckLumping", "%c lumping %d is not a divisor of %d, set to %d",
 257           which, lumping, nBins, l);
 258   return l;
 259 }
 260
 261 //____________________________________________________________________
 262 void
 263 AliPoissonCalculator::Reset(const TH2D* base)
 264 {
 265   //
 266   // Reset histogram
 267   //
 268   if (fBasic && fTotal && fEmpty) {
 269     fBasic->Reset();
 270     fTotal->Reset();
 271     fEmpty->Reset();
 272     return;
 273   }
 274
 275   if (!base) return;
 276
 277   Int_t    nXF   = base->GetNbinsX();
 278   Double_t xMin  = base->GetXaxis()->GetXmin();
 279   Double_t xMax  = base->GetXaxis()->GetXmax();
 280   Int_t    nYF   = base->GetNbinsY();
 281   Double_t yMin  = base->GetYaxis()->GetXmin();
 282   Double_t yMax  = base->GetYaxis()->GetXmax();
 283
 284   fXLumping = CheckLumping('X', nXF, fXLumping);
 285   fYLumping = CheckLumping('Y', nYF, fYLumping);
 286
 287   Int_t    nY    = nYF / fYLumping;
 288   Int_t    nX    = nXF / fXLumping;
 289
 290   if (fBasic != base) {
 291     fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
 292     fBasic->SetTitle(kBasicT);
 293     fBasic->SetDirectory(0);
 294     fBasic->Sumw2();
 295   }
 296
 297   fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax);
 298   fTotal->SetDirectory(0);
 299   fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle());
 300   fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle());
 301   fTotal->Sumw2();
 302
 303   fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
 304   fEmpty->SetTitle(kEmptyT);
 305   fEmpty->SetDirectory(0);
 306   // fEmpty->Sumw2();
 307 }
 308
 309 //____________________________________________________________________
 310 void
 311 AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight)
 312 {
 313   //
 314   // Fill in an observation
 315   //
 316   // Parameters:
 317   //    x       X value
 318   //    Y       Y value
 319   //    hit     True if hit
 320   //    weight  Weight if this
 321   //
 322   fTotal->Fill(x, y);
 323   if (hit) fBasic->Fill(x, y, weight);
 324   else     fEmpty->Fill(x, y);
 325 }
 326
 327 //____________________________________________________________________
 328 Double_t
 329 AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
 330 {
 331   if (total <= 0) return 0;
 332   if (empty < .001) empty = .001;
 333   return -TMath::Log(empty/total);
 334 }
 335 //____________________________________________________________________
 336 Double_t
 337 AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
 338 {
 339   if (total <= 0) return 0;
 340   if (TMath::Abs(empty-total) < .001) empty = total - .001;
 341   return 1 / (1 - empty / total);
 342 }
 343
 344 //____________________________________________________________________
 345 Int_t
 346 AliPoissonCalculator::GetReducedXBin(Int_t ix) const
 347 {
 348   if (!fBasic) return 0;
 349   Double_t mx = fBasic->GetXaxis()->GetBinCenter(ix);
 350   return GetReducedXBin(mx);
 351 }
 352 //____________________________________________________________________
 353 Int_t
 354 AliPoissonCalculator::GetReducedXBin(Double_t x) const
 355 {
 356   if (!fEmpty) return 0;
 357   return fEmpty->GetXaxis()->FindBin(x);
 358 }
 359 //____________________________________________________________________
 360 Int_t
 361 AliPoissonCalculator::GetReducedYBin(Int_t iy) const
 362 {
 363   if (!fBasic) return 0;
 364   Double_t my = fBasic->GetYaxis()->GetBinCenter(iy);
 365   return GetReducedYBin(my);
 366 }
 367 //____________________________________________________________________
 368 Int_t
 369 AliPoissonCalculator::GetReducedYBin(Double_t y) const
 370 {
 371   if (!fEmpty) return 0;
 372   return fEmpty->GetYaxis()->FindBin(y);
 373 }
 374
 375
 376
 377 //____________________________________________________________________
 378 TH2D*
 379 AliPoissonCalculator::Result(Bool_t correct)
 380 {
 381   //
 382   // Calculate result and store in @a output
 383   //
 384   // Return:
 385   //    The result histogram (fBase overwritten)
 386   //
 387
 388   // Double_t total = fXLumping * fYLumping;
 389
 390   for (Int_t ix = 1; ix <= fBasic->GetNbinsX(); ix++) {
 391     // Double_t x        = fBasic->GetXaxis()->GetBinCenter(ix);
 392     Int_t    jx       = GetReducedXBin(ix); // fEmpty->GetXaxis()->FindBin(x);
 393     for (Int_t iy = 1; iy <= fBasic->GetNbinsY(); iy++) {
 394       // Double_t y        = fBasic->GetYaxis()->GetBinCenter(iy);
 395       Int_t    jy       = GetReducedYBin(iy); // fEmpty->GetYaxis()->FindBin(y);
 396       Double_t empty    = fEmpty->GetBinContent(jx, jy);
 397       Double_t total    = fTotal->GetBinContent(jx, jy);
 398       Double_t hits     = fBasic->GetBinContent(ix,iy);
 399       // Mean in region of interest
 400       Double_t poissonM = CalculateMean(empty, total);
 401       Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1);
 402
 403       Double_t poissonV = hits * poissonM * poissonC;
 404       Double_t poissonE = TMath::Sqrt(poissonV);
 405       if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
 406
 407       fBasic->SetBinContent(ix,iy,poissonV);
 408       fBasic->SetBinError(ix,iy,poissonE);
 409     }
 410   }
 411   for (Int_t ix = 1; ix <= fEmpty->GetNbinsX(); ix++) {
 412     for (Int_t iy = 1; iy <= fEmpty->GetNbinsY(); iy++) {
 413       Double_t empty    = fEmpty->GetBinContent(ix, iy);
 414       Double_t total    = fTotal->GetBinContent(ix, iy);
 415       Double_t mean     = CalculateMean(empty, total);
 416       Double_t corr     = CalculateCorrection(empty, total);
 417       fEmptyVsTotal->Fill(total, empty);
 418       fMean->Fill(mean);
 419       fOcc->Fill(100 * (1 - empty/total));
 420       //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
 421       fCorr->Fill(mean, corr);
 422     }
 423   }
 424   return fBasic;
 425 }
 426
 427 //____________________________________________________________________
 428 void
 429 AliPoissonCalculator::Print(const Option_t*) const
 430 {
 431   //
 432   // Print information
 433   //
 434   // Parameters:
 435   //    option Not used
 436   //
 437   char ind[gROOT->GetDirLevel()+3];
 438   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 439   ind[gROOT->GetDirLevel()] = '\0';
 440   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 441             << std::boolalpha
 442             << ind << " X lumping:              " << fXLumping << '\n'
 443             << ind << " Y lumping:              " << fYLumping << '\n'
 444             << std::noboolalpha << std::endl;
 445 }
 446 //____________________________________________________________________
 447 void
 448 AliPoissonCalculator::Browse(TBrowser* b)
 449 {
 450   //
 451   // Browse this object
 452   //
 453   // Parameters:
 454   //    b Object to browse
 455   //
 456   if (fTotal)        b->Add(fTotal);
 457   if (fEmpty)        b->Add(fEmpty);
 458   if (fBasic)        b->Add(fBasic);
 459   if (fEmptyVsTotal) b->Add(fEmptyVsTotal);
 460   if (fMean)         b->Add(fMean);
 461   if (fOcc)          b->Add(fOcc);
 462   if (fCorr)         b->Add(fCorr);
 463 }
 464 //
 465 // EOF
 466 //