[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / DrawELossPoisson.C


Double_t
DrawRingELossPoisson(TList* p, UShort_t d, Char_t r, 
		     Double_t xmin, Double_t xmax)
{
  if (!p) return;

  TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r)));
  if (!ring) { 
    Error("DrawELossPoisson", "List FMD%d%c not found in %s",d,r,p->GetName());
    return;
  }
  
  TH2* corr = static_cast<TH2D*>(ring->FindObject("elossVsPoisson"));
  if (!corr) { 
    Error("DrawRingELossPoisson", "Histogram esdEloss not found in FMD%d%c",
	  d, r);
    return;
  }
  TPad* pad = static_cast<TPad*>(gPad);
  pad->SetGridy();
  pad->SetGridx();
  pad->SetLogz();
  if (d == 3) { 
    pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99, 
		 pad->GetYlowNDC()+pad->GetHNDC());
    pad->SetRightMargin(0.15);
  }
  pad->SetFillColor(0);
  if (xmax < 0) xmax = corr->GetXaxis()->GetXmax();
  corr->GetXaxis()->SetRangeUser(xmin,xmax);
  corr->GetYaxis()->SetRangeUser(xmin,xmax);
  corr->SetTitle(Form("FMD%d%c",d,r));
  Info("", "Entries: %d, integral: %f", corr->GetEntries(), corr->Integral()); 
  corr->Draw("colz");

  // Calculate the linear regression 
  Double_t dx    = (xmax-xmin);
  Double_t rxy   = corr->GetCorrelationFactor();
  Double_t sx    = corr->GetRMS(1);
  Double_t sy    = corr->GetRMS(2);
  Double_t sx2   = sx*sx;
  Double_t sy2   = sy*sy;
  Double_t cxy   = corr->GetCovariance();
  Double_t mx    = corr->GetMean(1);
  Double_t my    = corr->GetMean(2);
  Double_t delta = 1;
#if 0
  Double_t beta  = rxy * sy / sx;
#else
  Double_t beta  = ((sy2 - delta*sx2 + 
		     TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) + 
				 4*delta*cxy*cxy))
		    / 2 / cxy);
#endif
  Double_t alpha = my - beta * mx;
  TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax);
  f->SetParameters(alpha, beta);
  f->SetLineWidth(1);
  f->SetLineStyle(1);
  f->Draw("same");

  TLine* l = new TLine(xmin,xmin,xmax,xmax);
  l->SetLineWidth(1);
  l->SetLineStyle(2);
  l->SetLineColor(kBlack);
  l->Draw();
  // corr->GetXaxis()->SetRangeUser(0, 2);

  Info("", "FMD%d%c correlation coefficient: %9.5f%% "
       "line y = %f + %f * x", d, r, 100*rxy, alpha, beta);

  Double_t x = pad->GetLeftMargin()+.01;
  Double_t y = 1-pad->GetTopMargin()-.01;
  TLatex* ltx = new TLatex(x, y, Form("FMD%d%c", d, r));
  ltx->SetNDC();
  ltx->SetTextAlign(13);
  ltx->SetTextSize(0.08);
  ltx->Draw();
  y -= 0.12;
  ltx->SetTextSize(0.06);
  ltx->DrawLatex(x,y,"Deming regression: y=#alpha+#beta x");
  x += .02;
  y -= .06;
  ltx->SetTextSize(0.05);
  ltx->DrawLatex(x, y, Form("#alpha=%5.3f", alpha));
  y -= .06;
  ltx->DrawLatex(x, y, Form("#beta= %5.3f", beta));

  TLinearFitter* fitter = new TLinearFitter(1);
  fitter->SetFormula("1 ++ x");
  for (Int_t i = 1; i <= corr->GetNbinsX(); i++) { 
    x = corr->GetXaxis()->GetBinCenter(i);
    if (x < xmin || x > xmax) continue;
    for (Int_t j = 1; j <= corr->GetNbinsY(); j++) {
      y = corr->GetYaxis()->GetBinCenter(j);
      if (y < xmin || y > xmax) continue;
      Double_t c = corr->GetBinContent(i,j);
      if (c < .1) continue;
      fitter->AddPoint(&x, y, c);
    }
  }
  Bool_t robust = false;
  if (robust)
    fitter->EvalRobust();
  else 
    fitter->Eval();
  for (Int_t i = 0; i < 2; i++) { 
    std::cout << i << "\t" 
	      << fitter->GetParName(i) << "\t"
	      << fitter->GetParameter(i) << "\t";
    if (!robust) 
      std::cout << fitter->GetParError(i); 
    std::cout << std::endl;
  }
  Double_t chi2 = fitter->GetChisquare();
  Int_t    ndf  = (fitter->GetNpoints() - 
		   fitter->GetNumberFreeParameters() );
  std::cout << "chi2/ndf: " << chi2 << '/' << ndf 
	    << '=' << chi2 / ndf << std::endl;

  // corr->Scale(1. / corr->GetMaximum());
  pad->cd();
  return corr->GetCorrelationFactor();
}


void
DrawELossPoisson(const char* filename="forward.root", 
		 Double_t xmax=-1,
		 Double_t xmin=-1)
{
  gStyle->SetPalette(1);
  gStyle->SetOptFit(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetTitleW(.4);
  gStyle->SetTitleH(.1);
  gStyle->SetTitleX(.4);
  // gStyle->SetTitleY(.1);
  gStyle->SetTitleColor(0);
  gStyle->SetTitleStyle(0);
  gStyle->SetTitleBorderSize(0);
  
  TFile* file = TFile::Open(filename, "READ");
  if (!file) { 
    Error("DrawELossPoisson", "failed to open %s", filename);
    return;
  }

  TList* forward = static_cast<TList*>(file->Get("Forward"));
  if (!forward) { 
    Error("DrawELossPoisson", "List Forward not found in %s", filename);
    return;
  }

  TList* dc = static_cast<TList*>(forward->FindObject("fmdDensityCalculator"));
  if (!dc) { 
    Error("DrawELossPoisson", "List fmdDensityCalculator not found in Forward");
    return;
  }
  
  TCanvas* c = new TCanvas("elossVsPoisson", 
			   "N_ch from ELoss versus from Poisson", 900, 700);
  c->SetFillColor(0);
  c->SetBorderSize(0);
  c->SetBorderMode(0);
  c->SetHighLightColor(0);
  c->Divide(3, 2, 0, 0);

  Double_t corrs[5];
  c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax);
  c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax);
  c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax);
  c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax);
  c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax);

  TVirtualPad* p = c->cd(4);
  p->SetTopMargin(0.05);
  p->SetRightMargin(0.10);
  p->SetLeftMargin(0.15);
  p->SetBottomMargin(0.15);
  p->SetFillColor(0);

  TH1D* hc = new TH1D("corrs", "Correlation factors", 5, .5, 5.5);
  hc->SetFillColor(kRed+1);
  hc->SetFillStyle(3001);
  hc->SetMinimum(0.0);
  hc->SetMaximum(1.3);
  hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]);
  hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]);
  hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]);
  hc->GetXaxis()->SetBinLabel(4,"FMD3i"); hc->SetBinContent(4,corrs[3]);
  hc->GetXaxis()->SetBinLabel(5,"FMD3o"); hc->SetBinContent(5,corrs[4]);
  hc->GetXaxis()->SetLabelSize(0.08);
  hc->GetYaxis()->SetTitle("r");
  hc->SetMarkerSize(1.5);
  hc->Draw("text hist");

  // TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts"));
  // if (highCuts) highCuts->Draw("colz");
  c->cd();
  c->SaveAs("elossVsPoisson.png");
}
Commit	Line	Data
797161e8	1
0919aa56	2	Double_t
797161e8	3	DrawRingELossPoisson(TList* p, UShort_t d, Char_t r,
797161e8	4	Double_t xmin, Double_t xmax)
0919aa56	5	{
	6	if (!p) return;
	7
	8	TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r)));
	9	if (!ring) {
	10	Error("DrawELossPoisson", "List FMD%d%c not found in %s",d,r,p->GetName());
	11	return;
	12	}
	13
	14	TH2* corr = static_cast<TH2D*>(ring->FindObject("elossVsPoisson"));
	15	if (!corr) {
	16	Error("DrawRingELossPoisson", "Histogram esdEloss not found in FMD%d%c",
	17	d, r);
	18	return;
	19	}
b5c9a732	20	TPad* pad = static_cast<TPad*>(gPad);
	21	pad->SetGridy();
	22	pad->SetGridx();
	23	pad->SetLogz();
	24	if (d == 3) {
	25	pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99,
	26	pad->GetYlowNDC()+pad->GetHNDC());
	27	pad->SetRightMargin(0.15);
	28	}
	29	pad->SetFillColor(0);
797161e8	30	if (xmax < 0) xmax = corr->GetXaxis()->GetXmax();
	31	corr->GetXaxis()->SetRangeUser(xmin,xmax);
	32	corr->GetYaxis()->SetRangeUser(xmin,xmax);
0919aa56	33	corr->SetTitle(Form("FMD%d%c",d,r));
b5c9a732	34	Info("", "Entries: %d, integral: %f", corr->GetEntries(), corr->Integral());
0919aa56	35	corr->Draw("colz");
0919aa56	36
797161e8	37	// Calculate the linear regression
	38	Double_t dx = (xmax-xmin);
	39	Double_t rxy = corr->GetCorrelationFactor();
	40	Double_t sx = corr->GetRMS(1);
	41	Double_t sy = corr->GetRMS(2);
	42	Double_t sx2 = sx*sx;
	43	Double_t sy2 = sy*sy;
	44	Double_t cxy = corr->GetCovariance();
	45	Double_t mx = corr->GetMean(1);
	46	Double_t my = corr->GetMean(2);
	47	Double_t delta = 1;
	48	#if 0
	49	Double_t beta = rxy * sy / sx;
	50	#else
	51	Double_t beta = ((sy2 - delta*sx2 +
	52	TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) +
	53	4deltacxy*cxy))
	54	/ 2 / cxy);
	55	#endif
	56	Double_t alpha = my - beta * mx;
	57	TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax);
	58	f->SetParameters(alpha, beta);
	59	f->SetLineWidth(1);
	60	f->SetLineStyle(1);
	61	f->Draw("same");
	62
	63	TLine* l = new TLine(xmin,xmin,xmax,xmax);
0919aa56	64	l->SetLineWidth(1);
797161e8	65	l->SetLineStyle(2);
0919aa56	66	l->SetLineColor(kBlack);
	67	l->Draw();
	68	// corr->GetXaxis()->SetRangeUser(0, 2);
	69
797161e8	70	Info("", "FMD%d%c correlation coefficient: %9.5f%% "
	71	"line y = %f + %f * x", d, r, 100*rxy, alpha, beta);
	72
b5c9a732	73	Double_t x = pad->GetLeftMargin()+.01;
	74	Double_t y = 1-pad->GetTopMargin()-.01;
	75	TLatex* ltx = new TLatex(x, y, Form("FMD%d%c", d, r));
	76	ltx->SetNDC();
797161e8	77	ltx->SetTextAlign(13);
b5c9a732	78	ltx->SetTextSize(0.08);
797161e8	79	ltx->Draw();
b5c9a732	80	y -= 0.12;
	81	ltx->SetTextSize(0.06);
	82	ltx->DrawLatex(x,y,"Deming regression: y=#alpha+#beta x");
	83	x += .02;
	84	y -= .06;
797161e8	85	ltx->SetTextSize(0.05);
797161e8	86	ltx->DrawLatex(x, y, Form("#alpha=%5.3f", alpha));
b5c9a732	87	y -= .06;
797161e8	88	ltx->DrawLatex(x, y, Form("#beta= %5.3f", beta));
	89
	90	TLinearFitter* fitter = new TLinearFitter(1);
	91	fitter->SetFormula("1 ++ x");
	92	for (Int_t i = 1; i <= corr->GetNbinsX(); i++) {
	93	x = corr->GetXaxis()->GetBinCenter(i);
	94	if (x < xmin \|\| x > xmax) continue;
	95	for (Int_t j = 1; j <= corr->GetNbinsY(); j++) {
	96	y = corr->GetYaxis()->GetBinCenter(j);
	97	if (y < xmin \|\| y > xmax) continue;
	98	Double_t c = corr->GetBinContent(i,j);
	99	if (c < .1) continue;
	100	fitter->AddPoint(&x, y, c);
	101	}
	102	}
	103	Bool_t robust = false;
	104	if (robust)
	105	fitter->EvalRobust();
	106	else
	107	fitter->Eval();
	108	for (Int_t i = 0; i < 2; i++) {
	109	std::cout << i << "\t"
	110	<< fitter->GetParName(i) << "\t"
	111	<< fitter->GetParameter(i) << "\t";
	112	if (!robust)
	113	std::cout << fitter->GetParError(i);
	114	std::cout << std::endl;
	115	}
	116	Double_t chi2 = fitter->GetChisquare();
	117	Int_t ndf = (fitter->GetNpoints() -
	118	fitter->GetNumberFreeParameters() );
b5c9a732	119	std::cout << "chi2/ndf: " << chi2 << '/' << ndf
b5c9a732	120	<< '=' << chi2 / ndf << std::endl;
797161e8	121
b5c9a732	122	// corr->Scale(1. / corr->GetMaximum());
b5c9a732	123	pad->cd();
0919aa56	124	return corr->GetCorrelationFactor();
	125	}
	126
	127
	128	void
797161e8	129	DrawELossPoisson(const char* filename="forward.root",
797161e8	130	Double_t xmax=-1,
b5c9a732	131	Double_t xmin=-1)
0919aa56	132	{
	133	gStyle->SetPalette(1);
	134	gStyle->SetOptFit(0);
	135	gStyle->SetOptStat(0);
b5c9a732	136	gStyle->SetOptTitle(0);
0919aa56	137	gStyle->SetTitleW(.4);
0919aa56	138	gStyle->SetTitleH(.1);
b5c9a732	139	gStyle->SetTitleX(.4);
b5c9a732	140	// gStyle->SetTitleY(.1);
0919aa56	141	gStyle->SetTitleColor(0);
	142	gStyle->SetTitleStyle(0);
	143	gStyle->SetTitleBorderSize(0);
0919aa56	144
	145	TFile* file = TFile::Open(filename, "READ");
	146	if (!file) {
	147	Error("DrawELossPoisson", "failed to open %s", filename);
	148	return;
	149	}
	150
	151	TList* forward = static_cast<TList*>(file->Get("Forward"));
	152	if (!forward) {
	153	Error("DrawELossPoisson", "List Forward not found in %s", filename);
	154	return;
	155	}
	156
	157	TList* dc = static_cast<TList*>(forward->FindObject("fmdDensityCalculator"));
	158	if (!dc) {
	159	Error("DrawELossPoisson", "List fmdDensityCalculator not found in Forward");
	160	return;
	161	}
	162
	163	TCanvas* c = new TCanvas("elossVsPoisson",
	164	"N_ch from ELoss versus from Poisson", 900, 700);
	165	c->SetFillColor(0);
	166	c->SetBorderSize(0);
	167	c->SetBorderMode(0);
	168	c->SetHighLightColor(0);
	169	c->Divide(3, 2, 0, 0);
	170
	171	Double_t corrs[5];
797161e8	172	c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax);
	173	c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax);
	174	c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax);
	175	c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax);
	176	c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax);
0919aa56	177
	178	TVirtualPad* p = c->cd(4);
	179	p->SetTopMargin(0.05);
	180	p->SetRightMargin(0.10);
	181	p->SetLeftMargin(0.15);
	182	p->SetBottomMargin(0.15);
	183	p->SetFillColor(0);
	184
	185	TH1D* hc = new TH1D("corrs", "Correlation factors", 5, .5, 5.5);
	186	hc->SetFillColor(kRed+1);
	187	hc->SetFillStyle(3001);
	188	hc->SetMinimum(0.0);
797161e8	189	hc->SetMaximum(1.3);
0919aa56	190	hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]);
	191	hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]);
	192	hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]);
	193	hc->GetXaxis()->SetBinLabel(4,"FMD3i"); hc->SetBinContent(4,corrs[3]);
	194	hc->GetXaxis()->SetBinLabel(5,"FMD3o"); hc->SetBinContent(5,corrs[4]);
	195	hc->GetXaxis()->SetLabelSize(0.08);
	196	hc->GetYaxis()->SetTitle("r");
	197	hc->SetMarkerSize(1.5);
	198	hc->Draw("text hist");
	199
	200	// TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts"));
	201	// if (highCuts) highCuts->Draw("colz");
	202	c->cd();
797161e8	203	c->SaveAs("elossVsPoisson.png");
0919aa56	204	}
	205
	206
	207
	208