[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / tests / TestPoisson.C

/** 
 * 
 * 
 * @param nBins 
 * @param min 
 * @param max 
 *
 * @ingroup pwglf_forward_scripts_tests
 */
void
MakeIntegerAxis(Int_t& nBins, Double_t& min, Double_t& max)
{
  if (nBins <= 0) return;
  if (max <= 0)   max = nBins;
  
  Double_t dBin = max / nBins;
  min   = - dBin / 2;
  max   = max + dBin / 2;
  nBins = nBins + 1;
}

/** 
 * 
 * 
 * @param o 
 * @param useWeights 
 *
 * @ingroup pwglf_forward_scripts_tests
 */
void
TestPoisson(Double_t o=.3, bool useWeights=false, bool correct=true)
{
  const char* load = "$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C";
  gROOT->Macro(load);
  gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
  
  // --- Parameters of this script -----------------------------------
  Int_t      nBin =  5;  // Our detector matrix size 
  Int_t      nMax = TMath::Max(Int_t(nBin * nBin * o + .5)+nBin/2,nBin);  
  Int_t      nEv  = 10000; // Number of events
  Double_t   mp   = o;   // The 'hit' probability 


  TH2D* base = new TH2D("base", "Basic histogram", 
			nBin,-.5, nBin-.5, nBin, -.5, nBin-.5);
  base->SetXTitle("#eta");
  base->SetYTitle("#varphi");
  base->SetDirectory(0);
  base->SetOption("colz");

  Int_t tN1=nMax;    Double_t tMin1; Double_t tMax1;
  Int_t tN2=nMax*10; Double_t tMin2; Double_t tMax2=nMax;
  MakeIntegerAxis(tN1, tMin1, tMax1);
  MakeIntegerAxis(tN2, tMin2, tMax2);
  TH2D* corr = new TH2D("comp", "Comparison", 
			tN1, tMin1, tMax1, tN2, tMin2, tMax2);
  corr->SetXTitle("Input");
  corr->SetYTitle("Poisson");
  corr->SetDirectory(0);
  corr->SetOption("colz");
  corr->SetStats(0);
  TLine* lcorr = new TLine(0, 0, tMax2, tMax2);

  Int_t mm = TMath::Max(Int_t(nBin * o + .5),nBin/2);
  tN2=mm*10; tMax2 = mm;
  MakeIntegerAxis(tN2, tMin2, tMax2);
  Info("", "Making mean w/nbins=%d,range=[%f,%f]", tN2, tMin2, tMax2);
  TH2D* mean = new TH2D("mean", "Mean comparison", 
			tN2, tMin2, tMax2, tN2, tMin2, tMax2);
  mean->SetXTitle("Input");
  mean->SetYTitle("Poisson");
  mean->SetDirectory(0);
  mean->SetOption("colz");
  mean->SetStats(0);
  TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);

  TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
  dist->SetXTitle("s");
  dist->SetYTitle("P(s)");
  dist->SetFillColor(kRed+1);
  dist->SetFillStyle(3001);
  dist->SetDirectory(0);

  TH1D* diff = new TH1D("diff", "P-T", 100, -25, 25);
  diff->SetXTitle("Difference");
  diff->SetFillColor(kRed+1);
  diff->SetFillStyle(3001);
  diff->SetYTitle("Prob");

  AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
  c->Init(nBin ,nBin);

  for (Int_t i = 0; i < nEv; i++) { 
    c->Reset(base);
    base->Reset();

    for (Int_t iEta = 0; iEta < nBin; iEta++) { 
      for (Int_t iPhi = 0; iPhi < nBin; iPhi++) { 
	// Throw a die 
	Int_t m = gRandom->Poisson(mp);
	dist->Fill(m);

	// Fill into our base histogram 
	base->Fill(iEta, iPhi, m);

	// Fill into poisson calculator 
	c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
      }
    }
    // Calculate the result 
    TH2D* res = c->Result(correct);
    
    // Now loop and compare 
    Double_t mBase = 0;
    Double_t mPois = 0;
    for (Int_t iEta = 0; iEta < nBin; iEta++) { 
      for (Int_t iPhi = 0; iPhi < nBin; iPhi++) { 
	Double_t p = res->GetBinContent(iEta, iPhi);
	Double_t t = base->GetBinContent(iEta, iPhi);

	mBase += t;
	mPois += p;
	corr->Fill(t, p);
	diff->Fill(p-t);
      }
    }
    Int_t nn = nBin * nBin;
    mean->Fill(mBase / nn, mPois / nn);
  }

  TCanvas* cc = new TCanvas("c", "c", 900, 900);
  cc->SetFillColor(0);
  cc->SetFillStyle(0);
  cc->SetBorderMode(0);
  cc->SetRightMargin(0.02);
  cc->SetTopMargin(0.02);
  cc->Divide(2,2);
  
  TVirtualPad* pp = cc->cd(1);
  pp->SetFillColor(0);
  pp->SetFillStyle(0);
  pp->SetBorderMode(0);
  pp->SetRightMargin(0.15);
  pp->SetTopMargin(0.02);
  pp->SetLogz();
  pp->SetGridx();
  pp->SetGridy();
  corr->Draw();
  lcorr->Draw();

  pp = cc->cd(2);
  pp->SetFillColor(0);
  pp->SetFillStyle(0);
  pp->SetBorderMode(0);
  pp->SetRightMargin(0.02);
  pp->SetTopMargin(0.02);
#if 0
  c->GetMean()->Draw();
#elif 1 
  pp->SetLogy();
  diff->Draw();
#elif 1
  c->GetOccupancy()->Draw();
#else
  pp->SetLogy();
  dist->SetStats(0);
  dist->Scale(1. / dist->Integral());
  dist->Draw();
  TH1D* m1 = c->GetMean();
  m1->Scale(1. / m1->Integral());
  m1->Draw("same");
  Double_t eI;
  Double_t ii = 100 * dist->Integral(2, 0);
  TLatex* ll = new TLatex(.97, .85, 
			  Form("Input #bar{m}: %5.3f", mp));
  ll->SetNDC();
  ll->SetTextFont(132);
  ll->SetTextAlign(31);
  ll->Draw();
  ll->DrawLatex(.97, .75, Form("Result #bar{m}: %5.3f", dist->GetMean()));
  ll->DrawLatex(.97, .65, Form("Occupancy: #int_{1}^{#infty}P(s)ds = %6.2f%%",
			       ii));
			 
#endif

  pp = cc->cd(3);
  pp->SetFillColor(0);
  pp->SetFillStyle(0);
  pp->SetBorderMode(0);
  pp->SetRightMargin(0.15);
  pp->SetTopMargin(0.02);
  pp->SetGridx();
  pp->SetGridy();
  c->GetCorrection()->Draw();

  pp = cc->cd(4);
  pp->SetFillColor(0);
  pp->SetFillStyle(0);
  pp->SetBorderMode(0);
  pp->SetRightMargin(0.15);
  pp->SetTopMargin(0.02);
  pp->SetLogz();
  pp->SetGridx();
  pp->SetGridy();
  mean->Draw();
  lmean->Draw();

  cc->cd();
}

//
// EOF
//
Commit	Line	Data
56199f2b	1	/**
	2	*
	3	*
	4	* @param nBins
	5	* @param min
	6	* @param max
	7	*
bd6f5206	8	* @ingroup pwglf_forward_scripts_tests
56199f2b	9	*/
e5f38055	10	void
	11	MakeIntegerAxis(Int_t& nBins, Double_t& min, Double_t& max)
	12	{
	13	if (nBins <= 0) return;
	14	if (max <= 0) max = nBins;
	15
	16	Double_t dBin = max / nBins;
	17	min = - dBin / 2;
	18	max = max + dBin / 2;
	19	nBins = nBins + 1;
	20	}
	21
56199f2b	22	/**
	23	*
	24	*
	25	* @param o
	26	* @param useWeights
	27	*
bd6f5206	28	* @ingroup pwglf_forward_scripts_tests
56199f2b	29	*/
e5f38055	30	void
e18cb8bd	31	TestPoisson(Double_t o=.3, bool useWeights=false, bool correct=true)
e5f38055	32	{
bd6f5206	33	const char* load = "$ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/LoadLibs.C";
e5f38055	34	gROOT->Macro(load);
	35	gROOT->GetInterpreter()->UnloadFile(gSystem->ExpandPathName(load));
	36
	37	// --- Parameters of this script -----------------------------------
	38	Int_t nBin = 5; // Our detector matrix size
	39	Int_t nMax = TMath::Max(Int_t(nBin * nBin * o + .5)+nBin/2,nBin);
	40	Int_t nEv = 10000; // Number of events
	41	Double_t mp = o; // The 'hit' probability
	42
	43
	44	TH2D* base = new TH2D("base", "Basic histogram",
	45	nBin,-.5, nBin-.5, nBin, -.5, nBin-.5);
	46	base->SetXTitle("#eta");
	47	base->SetYTitle("#varphi");
	48	base->SetDirectory(0);
	49	base->SetOption("colz");
	50
	51	Int_t tN1=nMax; Double_t tMin1; Double_t tMax1;
	52	Int_t tN2=nMax*10; Double_t tMin2; Double_t tMax2=nMax;
	53	MakeIntegerAxis(tN1, tMin1, tMax1);
	54	MakeIntegerAxis(tN2, tMin2, tMax2);
	55	TH2D* corr = new TH2D("comp", "Comparison",
	56	tN1, tMin1, tMax1, tN2, tMin2, tMax2);
	57	corr->SetXTitle("Input");
	58	corr->SetYTitle("Poisson");
	59	corr->SetDirectory(0);
	60	corr->SetOption("colz");
	61	corr->SetStats(0);
	62	TLine* lcorr = new TLine(0, 0, tMax2, tMax2);
	63
	64	Int_t mm = TMath::Max(Int_t(nBin * o + .5),nBin/2);
	65	tN2=mm*10; tMax2 = mm;
	66	MakeIntegerAxis(tN2, tMin2, tMax2);
	67	Info("", "Making mean w/nbins=%d,range=[%f,%f]", tN2, tMin2, tMax2);
	68	TH2D* mean = new TH2D("mean", "Mean comparison",
	69	tN2, tMin2, tMax2, tN2, tMin2, tMax2);
	70	mean->SetXTitle("Input");
	71	mean->SetYTitle("Poisson");
	72	mean->SetDirectory(0);
	73	mean->SetOption("colz");
	74	mean->SetStats(0);
	75	TLine* lmean = new TLine(tMin2, tMin2, tMax2, tMax2);
	76
	77	TH1D* dist = new TH1D("dist", "Distribution of hits", tN1, tMin1, tMax1);
	78	dist->SetXTitle("s");
	79	dist->SetYTitle("P(s)");
	80	dist->SetFillColor(kRed+1);
	81	dist->SetFillStyle(3001);
	82	dist->SetDirectory(0);
e18cb8bd	83
	84	TH1D* diff = new TH1D("diff", "P-T", 100, -25, 25);
	85	diff->SetXTitle("Difference");
	86	diff->SetFillColor(kRed+1);
	87	diff->SetFillStyle(3001);
	88	diff->SetYTitle("Prob");
	89
e5f38055	90	AliPoissonCalculator* c = new AliPoissonCalculator("ignored");
	91	c->Init(nBin ,nBin);
	92
	93	for (Int_t i = 0; i < nEv; i++) {
	94	c->Reset(base);
	95	base->Reset();
	96
	97	for (Int_t iEta = 0; iEta < nBin; iEta++) {
	98	for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
	99	// Throw a die
	100	Int_t m = gRandom->Poisson(mp);
	101	dist->Fill(m);
	102
	103	// Fill into our base histogram
	104	base->Fill(iEta, iPhi, m);
	105
	106	// Fill into poisson calculator
	107	c->Fill(iEta, iPhi, m > 0, (useWeights ? m : 1));
	108	}
	109	}
	110	// Calculate the result
e18cb8bd	111	TH2D* res = c->Result(correct);
e5f38055	112
	113	// Now loop and compare
	114	Double_t mBase = 0;
	115	Double_t mPois = 0;
	116	for (Int_t iEta = 0; iEta < nBin; iEta++) {
	117	for (Int_t iPhi = 0; iPhi < nBin; iPhi++) {
	118	Double_t p = res->GetBinContent(iEta, iPhi);
	119	Double_t t = base->GetBinContent(iEta, iPhi);
	120
	121	mBase += t;
	122	mPois += p;
	123	corr->Fill(t, p);
e18cb8bd	124	diff->Fill(p-t);
e5f38055	125	}
	126	}
	127	Int_t nn = nBin * nBin;
	128	mean->Fill(mBase / nn, mPois / nn);
	129	}
	130
	131	TCanvas* cc = new TCanvas("c", "c", 900, 900);
	132	cc->SetFillColor(0);
	133	cc->SetFillStyle(0);
	134	cc->SetBorderMode(0);
	135	cc->SetRightMargin(0.02);
	136	cc->SetTopMargin(0.02);
	137	cc->Divide(2,2);
	138
	139	TVirtualPad* pp = cc->cd(1);
	140	pp->SetFillColor(0);
	141	pp->SetFillStyle(0);
	142	pp->SetBorderMode(0);
	143	pp->SetRightMargin(0.15);
	144	pp->SetTopMargin(0.02);
	145	pp->SetLogz();
	146	pp->SetGridx();
	147	pp->SetGridy();
	148	corr->Draw();
	149	lcorr->Draw();
	150
	151	pp = cc->cd(2);
	152	pp->SetFillColor(0);
	153	pp->SetFillStyle(0);
	154	pp->SetBorderMode(0);
	155	pp->SetRightMargin(0.02);
	156	pp->SetTopMargin(0.02);
e18cb8bd	157	#if 0
e5f38055	158	c->GetMean()->Draw();
e18cb8bd	159	#elif 1
	160	pp->SetLogy();
	161	diff->Draw();
e5f38055	162	#elif 1
	163	c->GetOccupancy()->Draw();
	164	#else
	165	pp->SetLogy();
	166	dist->SetStats(0);
	167	dist->Scale(1. / dist->Integral());
	168	dist->Draw();
	169	TH1D* m1 = c->GetMean();
	170	m1->Scale(1. / m1->Integral());
	171	m1->Draw("same");
	172	Double_t eI;
	173	Double_t ii = 100 * dist->Integral(2, 0);
	174	TLatex* ll = new TLatex(.97, .85,
	175	Form("Input #bar{m}: %5.3f", mp));
	176	ll->SetNDC();
	177	ll->SetTextFont(132);
	178	ll->SetTextAlign(31);
	179	ll->Draw();
	180	ll->DrawLatex(.97, .75, Form("Result #bar{m}: %5.3f", dist->GetMean()));
	181	ll->DrawLatex(.97, .65, Form("Occupancy: #int_{1}^{#infty}P(s)ds = %6.2f%%",
	182	ii));
	183
	184	#endif
	185
	186	pp = cc->cd(3);
	187	pp->SetFillColor(0);
	188	pp->SetFillStyle(0);
	189	pp->SetBorderMode(0);
	190	pp->SetRightMargin(0.15);
	191	pp->SetTopMargin(0.02);
	192	pp->SetGridx();
	193	pp->SetGridy();
	194	c->GetCorrection()->Draw();
	195
	196	pp = cc->cd(4);
	197	pp->SetFillColor(0);
	198	pp->SetFillStyle(0);
	199	pp->SetBorderMode(0);
	200	pp->SetRightMargin(0.15);
	201	pp->SetTopMargin(0.02);
	202	pp->SetLogz();
	203	pp->SetGridx();
	204	pp->SetGridy();
	205	mean->Draw();
	206	lmean->Draw();
	207
	208	cc->cd();
	209	}
	210
	211	//
	212	// EOF
	213	//
	214
	215