[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliPoissonCalculator.cxx

#include "AliPoissonCalculator.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$ as described by
// the 2D histogram passwd 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
// empty cells is then calculate in each region.  The mean
// multiplicity over the region is then determined as
//
// @f[
// \lange m\rangle = -\log\left(\frac{e}{t}\right)
// @f]
// where @f$ e@f$ is the number of empty cells and @f$t@f$ is the
// total number of cells in the region.  A correction for counting
// statistics, is then applied 
// @f{eqnarray*}
//    c &=& \frac{1}{1 - \exp{-\lange m\rangle}}
//      &=& \frac{1}{1 - \frac{e}{t}}
// @f{eqnarray*}
// and the final number in each cell is then 
// @f[
//   h_i c \lange m\rangle 
// @f] 
// where @f$h_i@f$ is the number of hits in the cell @f$i@f$ 
// 
//

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(),
    fXLumping(32), 
    fYLumping(4), 
    fTotal(0), 
    fEmpty(0), 
    fBasic(0),
    fEmptyVsTotal(0),
    fMean(0), 
    fOcc(0),
    fCorr(0)
{
  //
  // CTOR
  // 
}

//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator(const char*)
  : TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
    fXLumping(32), 
    fYLumping(4), 
    fTotal(0), 
    fEmpty(0), 
    fBasic(0),
    fEmptyVsTotal(0),
    fMean(0), 
    fOcc(0),
    fCorr(0)
{
  //
  // CTOR
  //
}
//____________________________________________________________________
AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
  : TNamed(o),
    fXLumping(o.fXLumping),
    fYLumping(o.fYLumping),
    fTotal(0), 
    fEmpty(0),
    fBasic(0), 
    fEmptyVsTotal(0),
    fMean(0), 
    fOcc(0),
    fCorr(0)
{
  Init();
  Reset(o.fBasic);
}

//____________________________________________________________________
AliPoissonCalculator::~AliPoissonCalculator()
{
  // CleanUp();
}

//____________________________________________________________________
void
AliPoissonCalculator::CleanUp()
{
  if (fTotal)        { delete fTotal;        fTotal        = 0; }
  if (fEmpty)        { delete fEmpty;        fEmpty        = 0; }
  if (fBasic)        { delete fBasic;        fBasic        = 0; }
  if (fEmptyVsTotal) { delete fEmptyVsTotal; fEmptyVsTotal = 0; }
  if (fMean)         { delete fMean;         fMean         = 0; }
  if (fOcc)          { delete fOcc;          fOcc          = 0; } 
  if (fCorr)         { delete fCorr;         fCorr         = 0; }
}
//____________________________________________________________________
AliPoissonCalculator&
AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
{
  if (&o == this) return *this;
  TNamed::operator=(o);
  fXLumping = o.fXLumping;
  fYLumping = o.fYLumping;
  CleanUp();
  Init();
  Reset(o.fBasic);
  return *this;
}

//____________________________________________________________________
void
AliPoissonCalculator::Init(Int_t xLumping, Int_t yLumping)
{
  // 
  // Initialize 
  // 
  if (xLumping > 0) fXLumping = xLumping;
  if (yLumping > 0) fYLumping = yLumping;

  //Create diagnostics if void
  if (fEmptyVsTotal) return;
  
  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", 
			   n, -.5, n-.5, n, -.5, n-.5);
  fEmptyVsTotal->SetXTitle("Total # bins");
  fEmptyVsTotal->SetYTitle("# empty bins");
  fEmptyVsTotal->SetZTitle("Correlation");
  fEmptyVsTotal->SetOption("colz");
  fEmptyVsTotal->SetDirectory(0);

  n = (fXLumping + fYLumping);
  fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
		   10*n+1, -.1, n+.1);
  fMean->SetXTitle("#bar{N_{ch}}=log(empty/total)");
  fMean->SetYTitle("Events");
  fMean->SetFillColor(kRed+1);
  fMean->SetFillStyle(3001);
  fMean->SetLineColor(kBlack);
  fMean->SetDirectory(0);

  fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
		  101, -.5, 100.5);
  fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
  fOcc->SetYTitle("Events");
  fOcc->SetFillColor(kBlue+1);
  fOcc->SetFillStyle(3001);
  fOcc->SetLineColor(kBlack);
  fOcc->SetDirectory(0);

  fCorr = new TH2D("correction", "Correction as function of mean N_{ch}", 
		   10*n+1, -.1, n+.1, 100, 0, 10);
  fCorr->SetXTitle("#bar{N_{ch}}");
  fCorr->SetYTitle("Correction 1/(1-e^{#bar{N_{c}}})");
  fCorr->SetZTitle("Events");
  fCorr->SetOption("colz");
  fCorr->SetDirectory(0); 
}
//____________________________________________________________________
void
AliPoissonCalculator::Output(TList* d)
{
  if (!d) return;
  if (!fEmptyVsTotal) MakeOutput();
  d->Add(fEmptyVsTotal);
  d->Add(fMean);
  d->Add(fOcc);
  d->Add(fCorr);
}

//____________________________________________________________________
void 
AliPoissonCalculator::SetLumping(UShort_t nx, UShort_t ny) 
{ 
  if (nx == fXLumping && ny == fYLumping && 
      fEmptyVsTotal && fTotal) 
    // Check if we have something to do. 
    return;
  CleanUp(); 
  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(const TH2D* base)
{
  // 
  // Reset histogram 
  // 
  if (fBasic && fTotal && fEmpty) {
    fBasic->Reset();
    fTotal->Reset();
    fEmpty->Reset();
    return;
  }

  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;
  Int_t    nX    = nXF / fXLumping;

  if (fBasic != base) { 
    fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
    fBasic->SetTitle(kBasicT);
    fBasic->SetDirectory(0);
    fBasic->Sumw2();
  }

  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(kEmptyN));
  fEmpty->SetTitle(kEmptyT);
  fEmpty->SetDirectory(0);
  // fEmpty->Sumw2();
}

//____________________________________________________________________
void
AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight)
{
  // 
  // Fill in an observation 
  // 
  // Parameters:
  //    x       X value 
  //    Y       Y value
  //    hit     True if hit 
  //    weight  Weight if this 
  //
  fTotal->Fill(x, y);
  if (hit) fBasic->Fill(x, y, weight);
  else     fEmpty->Fill(x, y);
}

//____________________________________________________________________
Double_t 
AliPoissonCalculator::CalculateMean(Double_t empty, Double_t total) const
{
  if (total <= 0) return 0;
  if (empty < .001) empty = .001;
  return -TMath::Log(empty/total);
}
//____________________________________________________________________
Double_t 
AliPoissonCalculator::CalculateCorrection(Double_t empty, Double_t total) const
{
  if (total <= 0) return 0;
  if (TMath::Abs(empty-total) < .001) empty = total - .001;
  return 1 / (1 - empty / total);
}

//____________________________________________________________________
Int_t
AliPoissonCalculator::GetReducedXBin(Int_t ix) const
{
  if (!fBasic) return 0;
  Double_t mx = fBasic->GetXaxis()->GetBinCenter(ix);
  return GetReducedXBin(mx);
}
//____________________________________________________________________
Int_t
AliPoissonCalculator::GetReducedXBin(Double_t x) const
{
  if (!fEmpty) return 0;
  return fEmpty->GetXaxis()->FindBin(x);
}
//____________________________________________________________________
Int_t
AliPoissonCalculator::GetReducedYBin(Int_t iy) const
{
  if (!fBasic) return 0;
  Double_t my = fBasic->GetYaxis()->GetBinCenter(iy);
  return GetReducedYBin(my);
}
//____________________________________________________________________
Int_t
AliPoissonCalculator::GetReducedYBin(Double_t y) const
{
  if (!fEmpty) return 0;
  return fEmpty->GetYaxis()->FindBin(y);
}


//____________________________________________________________________
TH2D*
AliPoissonCalculator::Result(Bool_t correct)
{
  // 
  // Calculate result and store in @a output
  // 
  // Return:
  //    The result histogram (fBase overwritten)
  //
  
  // Double_t total = fXLumping * fYLumping;
  
  for (Int_t ix = 1; ix <= fBasic->GetNbinsX(); ix++) { 
    // Double_t x        = fBasic->GetXaxis()->GetBinCenter(ix);
    Int_t    jx       = GetReducedXBin(ix); // fEmpty->GetXaxis()->FindBin(x);
    for (Int_t iy = 1; iy <= fBasic->GetNbinsY(); iy++) { 
      // Double_t y        = fBasic->GetYaxis()->GetBinCenter(iy);
      Int_t    jy       = GetReducedYBin(iy); // fEmpty->GetYaxis()->FindBin(y);
      Double_t empty    = fEmpty->GetBinContent(jx, jy);
      Double_t total    = fTotal->GetBinContent(jx, jy);
      Double_t hits     = fBasic->GetBinContent(ix,iy);
      // Mean in region of interest 
      Double_t poissonM = CalculateMean(empty, total);
      Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1);
      
      Double_t poissonV = hits * poissonM * poissonC;
      Double_t poissonE = TMath::Sqrt(poissonV);
      if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
	  
      fBasic->SetBinContent(ix,iy,poissonV);
      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 total    = fTotal->GetBinContent(ix, iy);
      Double_t mean     = CalculateMean(empty, total);
      Double_t corr     = CalculateCorrection(empty, total);
      fEmptyVsTotal->Fill(total, empty);
      fMean->Fill(mean);
      if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total));
      //Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
      fCorr->Fill(mean, corr);
    }
  }
}
//____________________________________________________________________
void
AliPoissonCalculator::Print(const Option_t*) const
{
  // 
  // Print information 
  // 
  // Parameters:
  //    option Not used
  //
  char ind[gROOT->GetDirLevel()+3];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << ClassName() << ": " << GetName() << '\n'
	    << std::boolalpha 
	    << ind << " X lumping:              " << fXLumping << '\n'
	    << ind << " Y lumping:              " << fYLumping << '\n'
	    << std::noboolalpha << std::endl;
}
//____________________________________________________________________
void
AliPoissonCalculator::Browse(TBrowser* b)
{
  // 
  // Browse this object
  // 
  // Parameters:
  //    b Object to browse 
  //
  if (fTotal)        b->Add(fTotal);
  if (fEmpty)        b->Add(fEmpty);
  if (fBasic)        b->Add(fBasic);
  if (fEmptyVsTotal) b->Add(fEmptyVsTotal);
  if (fMean)         b->Add(fMean);
  if (fOcc)          b->Add(fOcc);
  if (fCorr)         b->Add(fCorr);
}
// 
// EOF
//
Commit	Line	Data
71b70904	1	#include "AliPoissonCalculator.h"
8449e3e0	2	// #include "AliForwardCorrectionManager.h"
71b70904	3	#include <TH2D.h>
	4	#include <TBrowser.h>
	5	#include <TROOT.h>
ce63a99e	6	#include <TMath.h>
ce63a99e	7	#include <TList.h>
71b70904	8	#include <iostream>
d23503ee	9	#include <TAxis.h>
71b70904	10
ce63a99e	11	//
e18cb8bd	12	// A class to calculate the multiplicity in @f$(x,y)@f$ bins
ce63a99e	13	// using Poisson statistics.
ce63a99e	14	//
d23503ee	15	// The input is assumed to be binned in @f$(x,y)@f$ as described by
d23503ee	16	// the 2D histogram passwd to the Reset member function.
ce63a99e	17	//
ce63a99e	18	// The data is grouped in to regions as defined by the parameters
e18cb8bd	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
ce63a99e	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
d23503ee	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
71b70904	53	//____________________________________________________________________
	54	AliPoissonCalculator::AliPoissonCalculator()
	55	: TNamed(),
e18cb8bd	56	fXLumping(32),
e18cb8bd	57	fYLumping(4),
71b70904	58	fTotal(0),
71b70904	59	fEmpty(0),
ce63a99e	60	fBasic(0),
	61	fEmptyVsTotal(0),
	62	fMean(0),
	63	fOcc(0),
e18cb8bd	64	fCorr(0)
ce63a99e	65	{
	66	//
	67	// CTOR
	68	//
	69	}
71b70904	70
71b70904	71	//____________________________________________________________________
d23503ee	72	AliPoissonCalculator::AliPoissonCalculator(const char*)
71b70904	73	: TNamed("poissonCalculator", "Calculate N_ch using Poisson stat"),
e18cb8bd	74	fXLumping(32),
e18cb8bd	75	fYLumping(4),
71b70904	76	fTotal(0),
71b70904	77	fEmpty(0),
ce63a99e	78	fBasic(0),
	79	fEmptyVsTotal(0),
	80	fMean(0),
	81	fOcc(0),
e18cb8bd	82	fCorr(0)
ce63a99e	83	{
	84	//
	85	// CTOR
21d778b1	86	//
ce63a99e	87	}
	88	//____________________________________________________________________
	89	AliPoissonCalculator::AliPoissonCalculator(const AliPoissonCalculator& o)
	90	: TNamed(o),
e18cb8bd	91	fXLumping(o.fXLumping),
e18cb8bd	92	fYLumping(o.fYLumping),
ce63a99e	93	fTotal(0),
	94	fEmpty(0),
	95	fBasic(0),
	96	fEmptyVsTotal(0),
	97	fMean(0),
	98	fOcc(0),
e18cb8bd	99	fCorr(0)
ce63a99e	100	{
ce63a99e	101	Init();
d23503ee	102	Reset(o.fBasic);
ce63a99e	103	}
	104
	105	//____________________________________________________________________
	106	AliPoissonCalculator::~AliPoissonCalculator()
	107	{
08683480	108	// CleanUp();
ce63a99e	109	}
	110
	111	//____________________________________________________________________
	112	void
	113	AliPoissonCalculator::CleanUp()
	114	{
e18cb8bd	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; }
ce63a99e	122	}
	123	//____________________________________________________________________
	124	AliPoissonCalculator&
	125	AliPoissonCalculator::operator=(const AliPoissonCalculator& o)
	126	{
d015ecfe	127	if (&o == this) return *this;
ce63a99e	128	TNamed::operator=(o);
e18cb8bd	129	fXLumping = o.fXLumping;
e18cb8bd	130	fYLumping = o.fYLumping;
ce63a99e	131	CleanUp();
21d778b1	132	Init();
d23503ee	133	Reset(o.fBasic);
ce63a99e	134	return *this;
	135	}
	136
	137	//____________________________________________________________________
	138	void
e18cb8bd	139	AliPoissonCalculator::Init(Int_t xLumping, Int_t yLumping)
ce63a99e	140	{
	141	//
	142	// Initialize
	143	//
e18cb8bd	144	if (xLumping > 0) fXLumping = xLumping;
e18cb8bd	145	if (yLumping > 0) fYLumping = yLumping;
21d778b1	146
21d778b1	147	//Create diagnostics if void
ce63a99e	148	if (fEmptyVsTotal) return;
ce63a99e	149
d23503ee	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();
79909b8b	167
d23503ee	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);
79909b8b	180	}
d23503ee	181	//____________________________________________________________________
8449e3e0	182	void AliPoissonCalculator::MakeOutput()
8449e3e0	183	{
e18cb8bd	184	Int_t n = fXLumping * fYLumping + 1;
ce63a99e	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
e18cb8bd	194	n = (fXLumping + fYLumping);
ce63a99e	195	fMean = new TH1D("poissonMean", "Mean N_{ch} as calculated by Poisson",
ce63a99e	196	10*n+1, -.1, n+.1);
e18cb8bd	197	fMean->SetXTitle("#bar{N_{ch}}=log(empty/total)");
ce63a99e	198	fMean->SetYTitle("Events");
	199	fMean->SetFillColor(kRed+1);
	200	fMean->SetFillStyle(3001);
e2213ed5	201	fMean->SetLineColor(kBlack);
ce63a99e	202	fMean->SetDirectory(0);
	203
	204	fOcc = new TH1D("occupancy", "Occupancy = #int_{1}^{#infty}dN P(N)",
5ca83fee	205	101, -.5, 100.5);
ce63a99e	206	fOcc->SetXTitle("#int_{1}^{#infty}dN P(N) [%]");
	207	fOcc->SetYTitle("Events");
	208	fOcc->SetFillColor(kBlue+1);
	209	fOcc->SetFillStyle(3001);
e2213ed5	210	fOcc->SetLineColor(kBlack);
ce63a99e	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");
8449e3e0	219	fCorr->SetDirectory(0);
ce63a99e	220	}
ce63a99e	221	//____________________________________________________________________
	222	void
	223	AliPoissonCalculator::Output(TList* d)
	224	{
	225	if (!d) return;
79909b8b	226	if (!fEmptyVsTotal) MakeOutput();
ce63a99e	227	d->Add(fEmptyVsTotal);
	228	d->Add(fMean);
	229	d->Add(fOcc);
	230	d->Add(fCorr);
	231	}
71b70904	232
e18cb8bd	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
d23503ee	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
71b70904	259	//____________________________________________________________________
71b70904	260	void
d23503ee	261	AliPoissonCalculator::Reset(const TH2D* base)
71b70904	262	{
ce63a99e	263	//
	264	// Reset histogram
	265	//
e18cb8bd	266	if (fBasic && fTotal && fEmpty) {
71b70904	267	fBasic->Reset();
e18cb8bd	268	fTotal->Reset();
e18cb8bd	269	fEmpty->Reset();
71b70904	270	return;
71b70904	271	}
d23503ee	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
e18cb8bd	285	Int_t nY = nYF / fYLumping;
d23503ee	286	Int_t nX = nXF / fXLumping;
63d9ab53	287
d23503ee	288	if (fBasic != base) {
	289	fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
	290	fBasic->SetTitle(kBasicT);
	291	fBasic->SetDirectory(0);
	292	fBasic->Sumw2();
	293	}
21d778b1	294
d23503ee	295	fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax);
e18cb8bd	296	fTotal->SetDirectory(0);
	297	fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle());
	298	fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle());
	299	fTotal->Sumw2();
	300
d23503ee	301	fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
d23503ee	302	fEmpty->SetTitle(kEmptyT);
e18cb8bd	303	fEmpty->SetDirectory(0);
e18cb8bd	304	// fEmpty->Sumw2();
71b70904	305	}
	306
	307	//____________________________________________________________________
	308	void
e18cb8bd	309	AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight)
71b70904	310	{
ce63a99e	311	//
	312	// Fill in an observation
	313	//
	314	// Parameters:
e18cb8bd	315	// x X value
e18cb8bd	316	// Y Y value
ce63a99e	317	// hit True if hit
	318	// weight Weight if this
	319	//
e18cb8bd	320	fTotal->Fill(x, y);
	321	if (hit) fBasic->Fill(x, y, weight);
	322	else fEmpty->Fill(x, y);
71b70904	323	}
	324
	325	//____________________________________________________________________
ce63a99e	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
e18cb8bd	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
ce63a99e	375	//____________________________________________________________________
ce63a99e	376	TH2D*
e18cb8bd	377	AliPoissonCalculator::Result(Bool_t correct)
71b70904	378	{
ce63a99e	379	//
	380	// Calculate result and store in @a output
	381	//
	382	// Return:
	383	// The result histogram (fBase overwritten)
	384	//
21d778b1	385
e18cb8bd	386	// Double_t total = fXLumping * fYLumping;
21d778b1	387
e18cb8bd	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);
71b70904	397	// Mean in region of interest
ce63a99e	398	Double_t poissonM = CalculateMean(empty, total);
e18cb8bd	399	Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1);
21d778b1	400
ce63a99e	401	Double_t poissonV = hits * poissonM * poissonC;
ce63a99e	402	Double_t poissonE = TMath::Sqrt(poissonV);
71b70904	403	if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
71b70904	404
e18cb8bd	405	fBasic->SetBinContent(ix,iy,poissonV);
e18cb8bd	406	fBasic->SetBinError(ix,iy,poissonE);
71b70904	407	}
71b70904	408	}
77f97e3f CHC	409	return fBasic;
	410	}
	411
	412	//____________________________________________________________________
	413	void
	414	AliPoissonCalculator::FillDiagnostics()
	415	{
e18cb8bd	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);
ce63a99e	420	Double_t mean = CalculateMean(empty, total);
	421	Double_t corr = CalculateCorrection(empty, total);
	422	fEmptyVsTotal->Fill(total, empty);
	423	fMean->Fill(mean);
9efd555b	424	if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total));
21d778b1	425	//Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
ce63a99e	426	fCorr->Fill(mean, corr);
	427	}
	428	}
71b70904	429	}
71b70904	430	//____________________________________________________________________
	431	void
	432	AliPoissonCalculator::Print(const Option_t*) const
	433	{
ce63a99e	434	//
	435	// Print information
	436	//
	437	// Parameters:
	438	// option Not used
	439	//
71b70904	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
e18cb8bd	445	<< ind << " X lumping: " << fXLumping << '\n'
e18cb8bd	446	<< ind << " Y lumping: " << fYLumping << '\n'
9a621151	447	<< std::noboolalpha << std::endl;
71b70904	448	}
	449	//____________________________________________________________________
	450	void
	451	AliPoissonCalculator::Browse(TBrowser* b)
	452	{
ce63a99e	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);
71b70904	466	}
	467	//
	468	// EOF
	469	//