#include "AliPoissonCalculator.h"
-#include "AliForwardCorrectionManager.h"
+// #include "AliForwardCorrectionManager.h"
#include <TH2D.h>
#include <TBrowser.h>
#include <TROOT.h>
#include <TMath.h>
#include <TList.h>
#include <iostream>
+#include <TAxis.h>
//
// A class to calculate the multiplicity in @f$(x,y)@f$ bins
// using Poisson statistics.
//
-// The input is assumed to be binned in @f$(x,y)@f$ with the number of channels
-// input to the Reset member function.
+// The input is assumed to be binned in @f$(x,y)@f$ as described by
+// the 2D histogram passed to the Reset member function.
//
// The data is grouped in to regions as defined by the parameters
// fXLumping and fYLumping. The total number of cells and number of
//
//
+namespace {
+ const char* kBasicN = "basic";
+ const char* kEmptyN = "empty";
+ const char* kTotalN = "total";
+ const char* kBasicT = "Basic number of hits";
+ const char* kEmptyT = "Empty number of bins/region";
+ const char* kTotalT = "Total number of bins/region";
+}
+
+
+
+
//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator()
: TNamed(),
}
//____________________________________________________________________
-AliPoissonCalculator::AliPoissonCalculator(const char*/*, UShort_t d, Char_t r*/)
+AliPoissonCalculator::AliPoissonCalculator(const char*)
: TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
fXLumping(32),
fYLumping(4),
fCorr(0)
{
Init();
-// Reset();
+ Reset(o.fBasic);
}
//____________________________________________________________________
fYLumping = o.fYLumping;
CleanUp();
Init();
- // Reset();
+ Reset(o.fBasic);
return *this;
}
//Create diagnostics if void
if (fEmptyVsTotal) return;
- MakeOutput();
-
-
+ MakeOutput();
}
- //____________________________________________________________________
-void AliPoissonCalculator::MakeOutput() {
+//____________________________________________________________________
+void
+AliPoissonCalculator::Define(const TAxis& xaxis, const TAxis& yaxis)
+{
+ //
+ // Initialize
+ //
+ const Double_t* xBins = xaxis.GetXbins()->GetArray();
+ const Double_t* yBins = yaxis.GetXbins()->GetArray();
+ Int_t nX = xaxis.GetNbins();
+ Int_t nY = yaxis.GetNbins();
+ Double_t lX = xaxis.GetXmin();
+ Double_t hX = xaxis.GetXmax();
+ Double_t lY = yaxis.GetXmin();
+ Double_t hY = yaxis.GetXmax();
+
+ if (xBins) {
+ if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, yBins);
+ else fBasic = new TH2D(kBasicN, kBasicT, nX, xBins, nY, lY, hY);
+ }
+ else {
+ if (yBins) fBasic = new TH2D(kBasicN, kBasicT, nX, lX, hX, nY, yBins);
+ else fBasic = new TH2D(kBasicN, kBasicT, nX, lX, hX, nY, lY, hY);
+ }
+ fBasic->SetXTitle(xaxis.GetTitle());
+ fBasic->SetYTitle(yaxis.GetTitle());
+ Reset(fBasic);
+}
+//____________________________________________________________________
+void AliPoissonCalculator::MakeOutput()
+{
Int_t n = fXLumping * fYLumping + 1;
fEmptyVsTotal = new TH2D("emptyVsTotal",
"# of empty # bins vs total # bins",
fMean->SetDirectory(0);
fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
- 1000, 0, 100);
+ 101, -.5, 100.5);
fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
fOcc->SetYTitle("Events");
fOcc->SetFillColor(kBlue+1);
fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})");
fCorr->SetZTitle("Events");
fCorr->SetOption("colz");
- fCorr->SetDirectory(0);
-
-
+ fCorr->SetDirectory(0);
}
//____________________________________________________________________
void
Init(nx, ny);
}
+//____________________________________________________________________
+Int_t
+AliPoissonCalculator::CheckLumping(char which, Int_t nBins, Int_t lumping) const
+{
+ if ((nBins % lumping) == 0) return lumping;
+ Int_t l = lumping;
+ do {
+ l--;
+ } while (l > 0 && ((nBins % l) != 0));
+ Warning("CheckLumping", "%c lumping %d is not a divisor of %d, set to %d",
+ which, lumping, nBins, l);
+ return l;
+}
+
//____________________________________________________________________
void
-AliPoissonCalculator::Reset(UShort_t xChannels, UShort_t yChannels)
+AliPoissonCalculator::Reset(const TH2D* base)
{
//
// Reset histogram
//
-
if (fBasic && fTotal && fEmpty) {
fBasic->Reset();
fTotal->Reset();
fEmpty->Reset();
return;
}
- Int_t nXF = xChannels;
- Int_t nX = nXF / fXLumping;
- Double_t xMin = -0.5;
- Double_t xMax = nXF - 0.5;
- Int_t nYF = yChannels;
+
+ if (!base) return;
+
+ Int_t nXF = base->GetNbinsX();
+ Double_t xMin = base->GetXaxis()->GetXmin();
+ Double_t xMax = base->GetXaxis()->GetXmax();
+ Int_t nYF = base->GetNbinsY();
+ Double_t yMin = base->GetYaxis()->GetXmin();
+ Double_t yMax = base->GetYaxis()->GetXmax();
+
+ fXLumping = CheckLumping('X', nXF, fXLumping);
+ fYLumping = CheckLumping('Y', nYF, fYLumping);
+
Int_t nY = nYF / fYLumping;
- Double_t yMin = -0.5;
- Double_t yMax = nYF - 0.5;
+ Int_t nX = nXF / fXLumping;
- fBasic = new TH2D("basic", "basic number of hits",
- nX, xMin, xMax, nY, yMin, yMax);
- fBasic->SetTitle("Basic number of hits");
- fBasic->SetDirectory(0);
- fBasic->Sumw2();
+ if (fBasic != base) {
+ fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
+ fBasic->SetTitle(kBasicT);
+ fBasic->SetDirectory(0);
+ fBasic->Sumw2();
+ }
- fTotal = new TH2D("total", "Total number of bins/region",
- nX, xMin, xMax, nY, yMin, yMax);
+ fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax);
fTotal->SetDirectory(0);
fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle());
fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle());
fTotal->Sumw2();
- fEmpty = static_cast<TH2D*>(fTotal->Clone("empty"));
- fEmpty->SetTitle("Empty number of bins/region");
+ fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
+ fEmpty->SetTitle(kEmptyT);
fEmpty->SetDirectory(0);
// fEmpty->Sumw2();
}
fBasic->SetBinError(ix,iy,poissonE);
}
}
+ return fBasic;
+}
+
+//____________________________________________________________________
+void
+AliPoissonCalculator::FillDiagnostics()
+{
for (Int_t ix = 1; ix <= fEmpty->GetNbinsX(); ix++) {
for (Int_t iy = 1; iy <= fEmpty->GetNbinsY(); iy++) {
Double_t empty = fEmpty->GetBinContent(ix, iy);
Double_t corr = CalculateCorrection(empty, total);
fEmptyVsTotal->Fill(total, empty);
fMean->Fill(mean);
- fOcc->Fill(100 * (1 - empty/total));
+ if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total));
//Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
fCorr->Fill(mean, corr);
}
}
- return fBasic;
}
-
//____________________________________________________________________
void
AliPoissonCalculator::Print(const Option_t*) const
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff