[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / qa / DrawELossPoisson.C

/**
 * @file   DrawELossPoisson.C
 * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
 * @date   Thu Jul  7 10:24:58 2011
 * 
 * @brief  A script to draw the Poisson vs Energy Loss correlation 
 * 
 *
 * @deprecated Use QATrender instead 
 * @ingroup pwg2_forward_scripts_qa
 * 
 */
#ifndef __CINT__
# include <TH1.h>
# include <TH2.h>
# include <TList.h>
# include <TFile.h>
# include <TString.h>
# include <TError.h>
# include <TPad.h>
# include <TCanvas.h>
# include <TMath.h>
# include <TF1.h>
# include <TLine.h>
# include <TLatex.h>
# include <TLinearFitter.h>
# include <TStyle.h>
#else
class TList;
#endif

/** 
 * Draw the poisson @f$N_{ch}@f$ estimate against the @f$\Delta@f$
 * @f$N_{ch}@f$ estimate and do a regression line between the two 
 * 
 * @param p            List 
 * @param d            Detector
 * @param r            Ring
 * @param xmin         Minimum
 * @param xmax         Maximum
 * 
 * @return The regression coefficient 
 *
 * @deprecated Use QATrender instead 
 * @ingroup pwg2_forward_scripts_qa
 */
Double_t
DrawRingELossPoisson(TList* p, UShort_t d, Char_t r, 
		     Double_t xmin, Double_t xmax)
{
  if (!p) return 0;

  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 0;
  }
  
  TH2* corr = static_cast<TH2D*>(ring->FindObject("elossVsPoisson"));
  if (!corr) { 
    Error("DrawRingELossPoisson", "Histogram esdEloss not found in FMD%d%c",
	  d, r);
    return 0;
  }
  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", 
  //      int(corr->GetEntries()), corr->Integral()); 
  corr->Draw("colz");
  if (corr->GetEntries() <= 0) return 0;

  // 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
  if (TMath::Abs(cxy) < 1e-6) return 0;
  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);

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

  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();
#if 0
  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;
#endif

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

/** 
 * Draw the correlation between the Poisson @f$N_{ch}@f$ estimate
 * and the @f$\Delta@f$ @f$N_{ch}@f$ estimate and do a regression
 * line between the two for each ring
 * 
 * @param filename File to read
 * @param xmax     Minimum X
 * @param xmin     Maximum X 
 *
 * @deprecated Use QATrender instead 
 * @ingroup pwg2_forward_scripts_qa
 */
void
DrawELossPoisson(const char* filename="forward.root", 
		 const char* folder="ForwardResults",
		 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(folder));
  if (!forward) { 
    Error("DrawELossPoisson", "List %s not found in %s", folder, 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");
}
//
// EOF
//
Commit	Line	Data
1d86b7ca	1	/**
	2	* @file DrawELossPoisson.C
	3	* @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
	4	* @date Thu Jul 7 10:24:58 2011
	5	*
	6	* @brief A script to draw the Poisson vs Energy Loss correlation
	7	*
56199f2b	8	*
2dbde04b	9	* @deprecated Use QATrender instead
f0ef1e71	10	* @ingroup pwg2_forward_scripts_qa
1d86b7ca	11	*
1d86b7ca	12	*/
d015ecfe	13	#ifndef __CINT__
	14	# include <TH1.h>
	15	# include <TH2.h>
	16	# include <TList.h>
	17	# include <TFile.h>
	18	# include <TString.h>
	19	# include <TError.h>
	20	# include <TPad.h>
	21	# include <TCanvas.h>
	22	# include <TMath.h>
	23	# include <TF1.h>
	24	# include <TLine.h>
	25	# include <TLatex.h>
	26	# include <TLinearFitter.h>
	27	# include <TStyle.h>
	28	#else
	29	class TList;
	30	#endif
1d86b7ca	31
	32	/**
	33	* Draw the poisson @f$N_{ch}@f$ estimate against the @f$\Delta@f$
	34	* @f$N_{ch}@f$ estimate and do a regression line between the two
	35	*
	36	* @param p List
	37	* @param d Detector
	38	* @param r Ring
	39	* @param xmin Minimum
	40	* @param xmax Maximum
	41	*
	42	* @return The regression coefficient
56199f2b	43	*
2dbde04b	44	* @deprecated Use QATrender instead
f0ef1e71	45	* @ingroup pwg2_forward_scripts_qa
1d86b7ca	46	*/
0919aa56	47	Double_t
797161e8	48	DrawRingELossPoisson(TList* p, UShort_t d, Char_t r,
797161e8	49	Double_t xmin, Double_t xmax)
0919aa56	50	{
d015ecfe	51	if (!p) return 0;
0919aa56	52
	53	TList* ring = static_cast<TList*>(p->FindObject(Form("FMD%d%c",d,r)));
	54	if (!ring) {
	55	Error("DrawELossPoisson", "List FMD%d%c not found in %s",d,r,p->GetName());
d015ecfe	56	return 0;
0919aa56	57	}
	58
	59	TH2* corr = static_cast<TH2D*>(ring->FindObject("elossVsPoisson"));
	60	if (!corr) {
	61	Error("DrawRingELossPoisson", "Histogram esdEloss not found in FMD%d%c",
	62	d, r);
d015ecfe	63	return 0;
0919aa56	64	}
b5c9a732	65	TPad* pad = static_cast<TPad*>(gPad);
	66	pad->SetGridy();
	67	pad->SetGridx();
	68	pad->SetLogz();
	69	if (d == 3) {
	70	pad->SetPad(pad->GetXlowNDC(), pad->GetYlowNDC(), .99,
	71	pad->GetYlowNDC()+pad->GetHNDC());
	72	pad->SetRightMargin(0.15);
	73	}
	74	pad->SetFillColor(0);
797161e8	75	if (xmax < 0) xmax = corr->GetXaxis()->GetXmax();
	76	corr->GetXaxis()->SetRangeUser(xmin,xmax);
	77	corr->GetYaxis()->SetRangeUser(xmin,xmax);
0919aa56	78	corr->SetTitle(Form("FMD%d%c",d,r));
d015ecfe	79	// Info("", "Entries: %d, integral: %f",
d015ecfe	80	// int(corr->GetEntries()), corr->Integral());
0919aa56	81	corr->Draw("colz");
d015ecfe	82	if (corr->GetEntries() <= 0) return 0;
0919aa56	83
797161e8	84	// Calculate the linear regression
d015ecfe	85	// Double_t dx = (xmax-xmin);
797161e8	86	Double_t rxy = corr->GetCorrelationFactor();
	87	Double_t sx = corr->GetRMS(1);
	88	Double_t sy = corr->GetRMS(2);
	89	Double_t sx2 = sx*sx;
	90	Double_t sy2 = sy*sy;
	91	Double_t cxy = corr->GetCovariance();
	92	Double_t mx = corr->GetMean(1);
	93	Double_t my = corr->GetMean(2);
	94	Double_t delta = 1;
	95	#if 0
	96	Double_t beta = rxy * sy / sx;
	97	#else
d015ecfe	98	if (TMath::Abs(cxy) < 1e-6) return 0;
797161e8	99	Double_t beta = ((sy2 - delta*sx2 +
	100	TMath::Sqrt(TMath::Power(sy2-delta*sx2,2) +
	101	4deltacxy*cxy))
	102	/ 2 / cxy);
	103	#endif
	104	Double_t alpha = my - beta * mx;
	105	TF1* f = new TF1("f", "[0]+[1]*x", xmin, xmax);
	106	f->SetParameters(alpha, beta);
	107	f->SetLineWidth(1);
	108	f->SetLineStyle(1);
	109	f->Draw("same");
	110
	111	TLine* l = new TLine(xmin,xmin,xmax,xmax);
0919aa56	112	l->SetLineWidth(1);
797161e8	113	l->SetLineStyle(2);
0919aa56	114	l->SetLineColor(kBlack);
	115	l->Draw();
	116	// corr->GetXaxis()->SetRangeUser(0, 2);
	117
d015ecfe	118	#if 0
797161e8	119	Info("", "FMD%d%c correlation coefficient: %9.5f%% "
797161e8	120	"line y = %f + %f * x", d, r, 100*rxy, alpha, beta);
d015ecfe	121	#endif
797161e8	122
b5c9a732	123	Double_t x = pad->GetLeftMargin()+.01;
	124	Double_t y = 1-pad->GetTopMargin()-.01;
	125	TLatex* ltx = new TLatex(x, y, Form("FMD%d%c", d, r));
	126	ltx->SetNDC();
797161e8	127	ltx->SetTextAlign(13);
b5c9a732	128	ltx->SetTextSize(0.08);
797161e8	129	ltx->Draw();
b5c9a732	130	y -= 0.12;
	131	ltx->SetTextSize(0.06);
	132	ltx->DrawLatex(x,y,"Deming regression: y=#alpha+#beta x");
	133	x += .02;
	134	y -= .06;
797161e8	135	ltx->SetTextSize(0.05);
797161e8	136	ltx->DrawLatex(x, y, Form("#alpha=%5.3f", alpha));
b5c9a732	137	y -= .06;
797161e8	138	ltx->DrawLatex(x, y, Form("#beta= %5.3f", beta));
	139
	140	TLinearFitter* fitter = new TLinearFitter(1);
	141	fitter->SetFormula("1 ++ x");
	142	for (Int_t i = 1; i <= corr->GetNbinsX(); i++) {
	143	x = corr->GetXaxis()->GetBinCenter(i);
	144	if (x < xmin \|\| x > xmax) continue;
	145	for (Int_t j = 1; j <= corr->GetNbinsY(); j++) {
	146	y = corr->GetYaxis()->GetBinCenter(j);
	147	if (y < xmin \|\| y > xmax) continue;
	148	Double_t c = corr->GetBinContent(i,j);
	149	if (c < .1) continue;
	150	fitter->AddPoint(&x, y, c);
	151	}
	152	}
	153	Bool_t robust = false;
	154	if (robust)
	155	fitter->EvalRobust();
	156	else
	157	fitter->Eval();
d015ecfe	158	#if 0
797161e8	159	for (Int_t i = 0; i < 2; i++) {
	160	std::cout << i << "\t"
	161	<< fitter->GetParName(i) << "\t"
	162	<< fitter->GetParameter(i) << "\t";
	163	if (!robust)
	164	std::cout << fitter->GetParError(i);
	165	std::cout << std::endl;
	166	}
	167	Double_t chi2 = fitter->GetChisquare();
	168	Int_t ndf = (fitter->GetNpoints() -
	169	fitter->GetNumberFreeParameters() );
b5c9a732	170	std::cout << "chi2/ndf: " << chi2 << '/' << ndf
b5c9a732	171	<< '=' << chi2 / ndf << std::endl;
d015ecfe	172	#endif
797161e8	173
b5c9a732	174	// corr->Scale(1. / corr->GetMaximum());
b5c9a732	175	pad->cd();
0919aa56	176	return corr->GetCorrelationFactor();
	177	}
	178
1d86b7ca	179	/**
	180	* Draw the correlation between the Poisson @f$N_{ch}@f$ estimate
	181	* and the @f$\Delta@f$ @f$N_{ch}@f$ estimate and do a regression
	182	* line between the two for each ring
	183	*
	184	* @param filename File to read
	185	* @param xmax Minimum X
	186	* @param xmin Maximum X
56199f2b	187	*
2dbde04b	188	* @deprecated Use QATrender instead
f0ef1e71	189	* @ingroup pwg2_forward_scripts_qa
1d86b7ca	190	*/
0919aa56	191	void
797161e8	192	DrawELossPoisson(const char* filename="forward.root",
d015ecfe	193	const char* folder="ForwardResults",
797161e8	194	Double_t xmax=-1,
b5c9a732	195	Double_t xmin=-1)
0919aa56	196	{
	197	gStyle->SetPalette(1);
	198	gStyle->SetOptFit(0);
	199	gStyle->SetOptStat(0);
b5c9a732	200	gStyle->SetOptTitle(0);
0919aa56	201	gStyle->SetTitleW(.4);
0919aa56	202	gStyle->SetTitleH(.1);
b5c9a732	203	gStyle->SetTitleX(.4);
b5c9a732	204	// gStyle->SetTitleY(.1);
0919aa56	205	gStyle->SetTitleColor(0);
	206	gStyle->SetTitleStyle(0);
	207	gStyle->SetTitleBorderSize(0);
0919aa56	208
	209	TFile* file = TFile::Open(filename, "READ");
	210	if (!file) {
	211	Error("DrawELossPoisson", "failed to open %s", filename);
	212	return;
	213	}
	214
d015ecfe	215	TList* forward = static_cast<TList*>(file->Get(folder));
0919aa56	216	if (!forward) {
d015ecfe	217	Error("DrawELossPoisson", "List %s not found in %s", folder, filename);
0919aa56	218	return;
	219	}
	220
	221	TList* dc = static_cast<TList*>(forward->FindObject("fmdDensityCalculator"));
	222	if (!dc) {
	223	Error("DrawELossPoisson", "List fmdDensityCalculator not found in Forward");
	224	return;
	225	}
	226
	227	TCanvas* c = new TCanvas("elossVsPoisson",
	228	"N_ch from ELoss versus from Poisson", 900, 700);
	229	c->SetFillColor(0);
	230	c->SetBorderSize(0);
	231	c->SetBorderMode(0);
	232	c->SetHighLightColor(0);
	233	c->Divide(3, 2, 0, 0);
	234
	235	Double_t corrs[5];
797161e8	236	c->cd(1); corrs[0] = DrawRingELossPoisson(dc, 1, 'I', xmin, xmax);
	237	c->cd(2); corrs[1] = DrawRingELossPoisson(dc, 2, 'I', xmin, xmax);
	238	c->cd(5); corrs[2] = DrawRingELossPoisson(dc, 2, 'O', xmin, xmax);
	239	c->cd(3); corrs[3] = DrawRingELossPoisson(dc, 3, 'I', xmin, xmax);
	240	c->cd(6); corrs[4] = DrawRingELossPoisson(dc, 3, 'O', xmin, xmax);
0919aa56	241
	242	TVirtualPad* p = c->cd(4);
	243	p->SetTopMargin(0.05);
	244	p->SetRightMargin(0.10);
	245	p->SetLeftMargin(0.15);
	246	p->SetBottomMargin(0.15);
	247	p->SetFillColor(0);
	248
	249	TH1D* hc = new TH1D("corrs", "Correlation factors", 5, .5, 5.5);
	250	hc->SetFillColor(kRed+1);
	251	hc->SetFillStyle(3001);
	252	hc->SetMinimum(0.0);
797161e8	253	hc->SetMaximum(1.3);
0919aa56	254	hc->GetXaxis()->SetBinLabel(1,"FMD1i"); hc->SetBinContent(1,corrs[0]);
	255	hc->GetXaxis()->SetBinLabel(2,"FMD2i"); hc->SetBinContent(2,corrs[1]);
	256	hc->GetXaxis()->SetBinLabel(3,"FMD2o"); hc->SetBinContent(3,corrs[2]);
	257	hc->GetXaxis()->SetBinLabel(4,"FMD3i"); hc->SetBinContent(4,corrs[3]);
	258	hc->GetXaxis()->SetBinLabel(5,"FMD3o"); hc->SetBinContent(5,corrs[4]);
	259	hc->GetXaxis()->SetLabelSize(0.08);
	260	hc->GetYaxis()->SetTitle("r");
	261	hc->SetMarkerSize(1.5);
	262	hc->Draw("text hist");
	263
	264	// TH2D* highCuts = static_cast<TH2D*>(dc->FindObject("highCuts"));
	265	// if (highCuts) highCuts->Draw("colz");
	266	c->cd();
797161e8	267	c->SaveAs("elossVsPoisson.png");
0919aa56	268	}
1d86b7ca	269	//
	270	// EOF
	271	//