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