[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);
    }
  }
  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);
    }
  }
  return fBasic;
}
  
//____________________________________________________________________
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"
79909b8b	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	//____________________________________________________________________
79909b8b	182	void AliPoissonCalculator::MakeOutput() {
79909b8b	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");
	219	fCorr->SetDirectory(0);
21d778b1	220
79909b8b	221
ce63a99e	222	}
ce63a99e	223	//____________________________________________________________________
	224	void
	225	AliPoissonCalculator::Output(TList* d)
	226	{
	227	if (!d) return;
79909b8b	228	if (!fEmptyVsTotal) MakeOutput();
ce63a99e	229	d->Add(fEmptyVsTotal);
	230	d->Add(fMean);
	231	d->Add(fOcc);
	232	d->Add(fCorr);
	233	}
71b70904	234
e18cb8bd	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
d23503ee	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
71b70904	261	//____________________________________________________________________
71b70904	262	void
d23503ee	263	AliPoissonCalculator::Reset(const TH2D* base)
71b70904	264	{
ce63a99e	265	//
	266	// Reset histogram
	267	//
e18cb8bd	268	if (fBasic && fTotal && fEmpty) {
71b70904	269	fBasic->Reset();
e18cb8bd	270	fTotal->Reset();
e18cb8bd	271	fEmpty->Reset();
71b70904	272	return;
71b70904	273	}
d23503ee	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
e18cb8bd	287	Int_t nY = nYF / fYLumping;
d23503ee	288	Int_t nX = nXF / fXLumping;
63d9ab53	289
d23503ee	290	if (fBasic != base) {
	291	fBasic = static_cast<TH2D*>(base->Clone(kBasicN));
	292	fBasic->SetTitle(kBasicT);
	293	fBasic->SetDirectory(0);
	294	fBasic->Sumw2();
	295	}
21d778b1	296
d23503ee	297	fTotal = new TH2D(kTotalN, kTotalT, nX, xMin, xMax, nY, yMin, yMax);
e18cb8bd	298	fTotal->SetDirectory(0);
	299	fTotal->SetXTitle(fBasic->GetXaxis()->GetTitle());
	300	fTotal->SetYTitle(fBasic->GetYaxis()->GetTitle());
	301	fTotal->Sumw2();
	302
d23503ee	303	fEmpty = static_cast<TH2D*>(fTotal->Clone(kEmptyN));
d23503ee	304	fEmpty->SetTitle(kEmptyT);
e18cb8bd	305	fEmpty->SetDirectory(0);
e18cb8bd	306	// fEmpty->Sumw2();
71b70904	307	}
	308
	309	//____________________________________________________________________
	310	void
e18cb8bd	311	AliPoissonCalculator::Fill(UShort_t x, UShort_t y, Bool_t hit, Double_t weight)
71b70904	312	{
ce63a99e	313	//
	314	// Fill in an observation
	315	//
	316	// Parameters:
e18cb8bd	317	// x X value
e18cb8bd	318	// Y Y value
ce63a99e	319	// hit True if hit
	320	// weight Weight if this
	321	//
e18cb8bd	322	fTotal->Fill(x, y);
	323	if (hit) fBasic->Fill(x, y, weight);
	324	else fEmpty->Fill(x, y);
71b70904	325	}
	326
	327	//____________________________________________________________________
ce63a99e	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
e18cb8bd	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
ce63a99e	377	//____________________________________________________________________
ce63a99e	378	TH2D*
e18cb8bd	379	AliPoissonCalculator::Result(Bool_t correct)
71b70904	380	{
ce63a99e	381	//
	382	// Calculate result and store in @a output
	383	//
	384	// Return:
	385	// The result histogram (fBase overwritten)
	386	//
21d778b1	387
e18cb8bd	388	// Double_t total = fXLumping * fYLumping;
21d778b1	389
e18cb8bd	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);
71b70904	399	// Mean in region of interest
ce63a99e	400	Double_t poissonM = CalculateMean(empty, total);
e18cb8bd	401	Double_t poissonC = (correct ? CalculateCorrection(empty, total) : 1);
21d778b1	402
ce63a99e	403	Double_t poissonV = hits * poissonM * poissonC;
ce63a99e	404	Double_t poissonE = TMath::Sqrt(poissonV);
71b70904	405	if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
71b70904	406
e18cb8bd	407	fBasic->SetBinContent(ix,iy,poissonV);
e18cb8bd	408	fBasic->SetBinError(ix,iy,poissonE);
71b70904	409	}
71b70904	410	}
e18cb8bd	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);
ce63a99e	415	Double_t mean = CalculateMean(empty, total);
	416	Double_t corr = CalculateCorrection(empty, total);
	417	fEmptyVsTotal->Fill(total, empty);
	418	fMean->Fill(mean);
9efd555b	419	if (total > 1e-6) fOcc->Fill(100 * (1 - empty/total));
21d778b1	420	//Old fOcc->Fill(100 * (1 - TMath::PoissonI(0,mean)));
ce63a99e	421	fCorr->Fill(mean, corr);
	422	}
	423	}
	424	return fBasic;
71b70904	425	}
	426
	427	//____________________________________________________________________
	428	void
	429	AliPoissonCalculator::Print(const Option_t*) const
	430	{
ce63a99e	431	//
	432	// Print information
	433	//
	434	// Parameters:
	435	// option Not used
	436	//
71b70904	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
e18cb8bd	442	<< ind << " X lumping: " << fXLumping << '\n'
e18cb8bd	443	<< ind << " Y lumping: " << fYLumping << '\n'
9a621151	444	<< std::noboolalpha << std::endl;
71b70904	445	}
	446	//____________________________________________________________________
	447	void
	448	AliPoissonCalculator::Browse(TBrowser* b)
	449	{
ce63a99e	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);
71b70904	463	}
	464	//
	465	// EOF
	466	//