[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / qa / DrawMCResult.C

THStack*
GetStack(const TList& forward,  const char* sub, const char* name)
{
  TList* lsub = static_cast<TList*>(forward.FindObject(sub));
  if (!lsub) { 
    Warning("GetStack", "Sub list %s not found in %s", sub, forward.GetName());
    return 0;
  }
  THStack* ret = static_cast<THStack*>(lsub->FindObject(name));
  if (!ret) 
    Warning("GetStack" "Stack %s not found in %s", name, sub);
  return ret;
}

TH1* 
Rebin(TH1* h, Int_t rebin)
{
  if (rebin <= 1) return h;
  h->Rebin(rebin);
  h->Scale(1. / rebin);
  return h;
}

TH1*
Ratio(const TH1* h1, const TH1* h2)
{
  if (!h1) return;
  if (!h2) return;
  
  TH1* copy = static_cast<TH1*>(h2->Clone("tmp"));
  copy->SetName(Form("%s_%s", h2->GetName(), h1->GetName()));
  copy->SetTitle(Form("%s/%s", h2->GetTitle(), h1->GetTitle()));
  copy->SetDirectory(0);
  copy->Divide(h1);

  return copy;
}

Int_t 
Ratio(THStack* r, const THStack* h1, const THStack* h2)
{
  if (!h1) return 0;
  if (!h2) return 0;

  int n1 = h1->GetHists()->GetEntries();
  int n2 = h2->GetHists()->GetEntries();
  int nH = 0;
  for (int i = 0; i < n1 && i < n2; i++) { 
    TH1* hh1 = static_cast<TH1*>(h1->GetHists()->At(i));
    TH1* hh2 = static_cast<TH1*>(h2->GetHists()->At(i));
    TH1* h   = Ratio(hh1, hh2);
    if (!h) continue;
    nH++;
    r->Add(h);
  }
  return nH;
}

void
DrawMCResult(const char* filename="forward.root", Int_t rebin=1,
	     Bool_t ratios=true)
{
  gStyle->SetPalette(1);
  gStyle->SetOptFit(0);
  gStyle->SetOptStat(0);

  TFile* file = TFile::Open(filename, "READ");
  if (!file) { 
    Error("DrawMCResult", "failed to open %s", filename);
    return;
  }

  TList* forward = static_cast<TList*>(file->Get("Forward"));
  if (!forward) { 
    Error("DrawMCResult", "List Forward not found in %s", filename);
    return;
  }
  THStack* res    = GetStack(*forward, "ringResults", "all");
  THStack* mcRes  = GetStack(*forward, "mcRingResults", "all");
  THStack* deltas = GetStack(*forward, "fmdSharingFilter", "sums");
  THStack* nchs   = GetStack(*forward, "fmdDensityCalculator", "sums");
  THStack* prims  = GetStack(*forward, "fmdCorrector", "sums");
  
  TH1* sumEta = static_cast<TH1*>(forward->FindObject("mcSumEta"));
  if (!sumEta) { 
    Warning("DrawMCResults", "mcSumEta not found in Forward");
  }


  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetTitleStyle(0);
  // gStyle->SetTitleColor(kBlack);


  TCanvas* c = new TCanvas("c", "C", 900, 700);
  c->SetFillColor(0);
  c->SetBorderSize(0);
  c->SetTopMargin(0.05);
  c->SetRightMargin(0.05);

  Double_t y1 = (ratios ? .3 : 0);
  TPad* p1 = new TPad("p1", "p1", 0, y1, 1, 1, 0, 0, 0); 
  p1->SetBottomMargin(ratios ? 0.01 : .10);
  p1->SetFillColor(0);
  p1->SetBorderSize(0);
  p1->SetTopMargin(0.05);
  p1->SetRightMargin(0.05);
  p1->Draw();
  p1->cd();

  THStack* all = new THStack("all", "Analysis steps");
  if (res) {
    res->SetTitle("dN_{ch}/d#eta");
    TIter next(res->GetHists());
    TH1* h = 0;
    while ((h = static_cast<TH1*>(next()))) all->Add(Rebin(h,rebin));
  }
  if (mcRes) {
    mcRes->SetTitle("Track-Refs");
    TIter next(mcRes->GetHists());
    TH1* h = 0;
    while ((h = static_cast<TH1*>(next()))) all->Add(Rebin(h,rebin));
  }
  if (deltas) {
    deltas->SetTitle("#sum_{} #Delta/#Delta_{mip}");
    TIter next(deltas->GetHists());
    TH1* h = 0;
    while ((h = static_cast<TH1*>(next()))) { 
      h->SetMarkerStyle(25);
      all->Add(Rebin(h,rebin));
    }
  }
  if (nchs) {
    nchs->SetTitle("#sum_{} N_{ch,incl}");
    TIter next(nchs->GetHists());
    TH1* h = 0;
    while ((h = static_cast<TH1*>(next()))) { 
      h->SetMarkerStyle(21);
      all->Add(Rebin(h,rebin));
    }
  }
  if (prims) {
    prims->SetTitle("#sum_{} N_{ch,primary}");
    TIter next(prims->GetHists());
    TH1* h = 0;
    while ((h = static_cast<TH1*>(next()))) { 
      h->SetMarkerStyle(22);
      all->Add(Rebin(h,rebin));
    }
  }
  if (sumEta) all->Add(sumEta);
  all->Draw("nostack");
  all->GetHistogram()->SetXTitle("#eta");
  all->GetHistogram()->SetYTitle("signal");
  all->GetHistogram()->GetXaxis()->SetLabelFont(132);
  all->GetHistogram()->GetXaxis()->SetTitleFont(132);
  all->GetHistogram()->GetYaxis()->SetLabelFont(132);
  all->GetHistogram()->GetYaxis()->SetTitleFont(132);
  c->SetGridx();

  TLegend* l = new TLegend(.33, .2, .53, .9);
  TLegendEntry* e = 0;
  l->SetFillColor(0);
  l->SetFillStyle(0);
  l->SetBorderSize(0);
  l->SetNColumns(1);
  l->SetTextFont(132);
  Int_t i = 0;
  if (res) { 
    TIter next(res->GetHists());
    TH1* h = 0;
    while ((h = static_cast<TH1*>(next()))) {
      e = l->AddEntry(Form("dummy%02d", i++),h->GetTitle(),"pl");
      e->SetMarkerStyle(20);
      e->SetMarkerColor(h->GetMarkerColor());
    }
    e = l->AddEntry(Form("dummy%02d", i++),res->GetTitle(),"pl");
    e->SetMarkerStyle(20);
    e->SetMarkerColor(kBlack);
  }
  if (deltas) { 
    e = l->AddEntry(Form("dummy%02d", i++), deltas->GetTitle(),"pl");
    TH1* h = static_cast<TH1*>(deltas->GetHists()->At(0));
    e->SetMarkerStyle(h->GetMarkerStyle());
    e->SetMarkerColor(kBlack);
  }
  if (nchs) { 
    e = l->AddEntry(Form("dummy%02d",i++),nchs->GetTitle(),"pl");
    TH1* h = static_cast<TH1*>(nchs->GetHists()->At(0));
    e->SetMarkerStyle(h->GetMarkerStyle());
    e->SetMarkerColor(kBlack);
  }
  if (prims) { 
    e = l->AddEntry(Form("dummy%02d", i++), prims->GetTitle(),"pl");
    TH1* h = static_cast<TH1*>(prims->GetHists()->At(0));
    e->SetMarkerStyle(h->GetMarkerStyle());
    e->SetMarkerColor(kBlack);
  }

  if (mcRes) { 
    e = l->AddEntry(Form("dummy%02d", i++), mcRes->GetTitle(), "pl");
    TH1* h = static_cast<TH1*>(mcRes->GetHists()->At(0));
    e->SetMarkerStyle(h->GetMarkerStyle());
    e->SetMarkerColor(kBlack);
  }

  if (sumEta) l->AddEntry(sumEta);
  l->Draw();


  if (!ratios) return;

  c->cd();
  TPad* p2 = new TPad("p2", "p2", 0, 0, 1, y1, 0, 0, 0); 
  p2->SetTopMargin(0);
  p2->SetFillColor(0);
  p2->SetBorderSize(0);
  p2->SetRightMargin(0.05);
  p2->Draw();
  p2->cd();

  THStack* rs = new THStack("ratios", "Ratios");
  Int_t nDN = Ratio(rs, deltas, nchs);
  Int_t nNR = Ratio(rs, nchs, res);
  Int_t nRP = Ratio(rs, res, prims);
  rs->Draw("nostack");

  TLegend* ll = new TLegend(.38, .2, .48, .9);
  ll->SetFillColor(0);
  ll->SetFillStyle(0);
  ll->SetBorderSize(0);
  ll->SetNColumns(1);
  ll->SetTextFont(132);
  if (nDN) {
    e = ll->AddEntry("d1", Form("#frac{%s}{%s}", 
					      nchs->GetTitle(), 
					      deltas->GetTitle()), "lp");
    e->SetMarkerStyle(21);
  }
  if (nNR) {
    e = ll->AddEntry("d2", Form("#frac{%s}{%s}", res->GetTitle(), 
				nchs->GetTitle()), "lp");
    e->SetMarkerStyle(20);
  }
  if (nRP) {
    e = ll->AddEntry("d3", Form("#frac{%s}{%s}", prims->GetTitle(), 
			    res->GetTitle()), "lp");
    e->SetMarkerStyle(22);
  }
  ll->Draw();
  

}
Commit	Line	Data
0919aa56	1	THStack*
	2	GetStack(const TList& forward, const char* sub, const char* name)
	3	{
	4	TList* lsub = static_cast<TList*>(forward.FindObject(sub));
	5	if (!lsub) {
	6	Warning("GetStack", "Sub list %s not found in %s", sub, forward.GetName());
	7	return 0;
	8	}
	9	THStack* ret = static_cast<THStack*>(lsub->FindObject(name));
	10	if (!ret)
	11	Warning("GetStack" "Stack %s not found in %s", name, sub);
	12	return ret;
	13	}
	14
	15	TH1*
	16	Rebin(TH1* h, Int_t rebin)
	17	{
	18	if (rebin <= 1) return h;
	19	h->Rebin(rebin);
	20	h->Scale(1. / rebin);
	21	return h;
	22	}
	23
8b1c306c	24	TH1*
	25	Ratio(const TH1* h1, const TH1* h2)
	26	{
	27	if (!h1) return;
	28	if (!h2) return;
	29
	30	TH1* copy = static_cast<TH1*>(h2->Clone("tmp"));
	31	copy->SetName(Form("%s_%s", h2->GetName(), h1->GetName()));
	32	copy->SetTitle(Form("%s/%s", h2->GetTitle(), h1->GetTitle()));
	33	copy->SetDirectory(0);
	34	copy->Divide(h1);
	35
	36	return copy;
	37	}
	38
	39	Int_t
	40	Ratio(THStack* r, const THStack* h1, const THStack* h2)
	41	{
	42	if (!h1) return 0;
	43	if (!h2) return 0;
	44
	45	int n1 = h1->GetHists()->GetEntries();
	46	int n2 = h2->GetHists()->GetEntries();
	47	int nH = 0;
	48	for (int i = 0; i < n1 && i < n2; i++) {
	49	TH1* hh1 = static_cast<TH1*>(h1->GetHists()->At(i));
	50	TH1* hh2 = static_cast<TH1*>(h2->GetHists()->At(i));
	51	TH1* h = Ratio(hh1, hh2);
	52	if (!h) continue;
	53	nH++;
	54	r->Add(h);
	55	}
	56	return nH;
	57	}
	58
0919aa56	59	void
8b1c306c	60	DrawMCResult(const char* filename="forward.root", Int_t rebin=1,
8b1c306c	61	Bool_t ratios=true)
0919aa56	62	{
	63	gStyle->SetPalette(1);
	64	gStyle->SetOptFit(0);
	65	gStyle->SetOptStat(0);
	66
	67	TFile* file = TFile::Open(filename, "READ");
	68	if (!file) {
	69	Error("DrawMCResult", "failed to open %s", filename);
	70	return;
	71	}
	72
	73	TList* forward = static_cast<TList*>(file->Get("Forward"));
	74	if (!forward) {
	75	Error("DrawMCResult", "List Forward not found in %s", filename);
	76	return;
	77	}
	78	THStack* res = GetStack(*forward, "ringResults", "all");
	79	THStack* mcRes = GetStack(*forward, "mcRingResults", "all");
	80	THStack* deltas = GetStack(*forward, "fmdSharingFilter", "sums");
	81	THStack* nchs = GetStack(*forward, "fmdDensityCalculator", "sums");
	82	THStack* prims = GetStack(*forward, "fmdCorrector", "sums");
	83
	84	TH1* sumEta = static_cast<TH1*>(forward->FindObject("mcSumEta"));
	85	if (!sumEta) {
	86	Warning("DrawMCResults", "mcSumEta not found in Forward");
	87	}
	88
	89
8b1c306c	90	gStyle->SetTitleBorderSize(0);
	91	gStyle->SetTitleFillColor(0);
	92	gStyle->SetTitleStyle(0);
	93	// gStyle->SetTitleColor(kBlack);
	94
0919aa56	95
	96	TCanvas* c = new TCanvas("c", "C", 900, 700);
	97	c->SetFillColor(0);
	98	c->SetBorderSize(0);
	99	c->SetTopMargin(0.05);
	100	c->SetRightMargin(0.05);
8b1c306c	101
	102	Double_t y1 = (ratios ? .3 : 0);
	103	TPad* p1 = new TPad("p1", "p1", 0, y1, 1, 1, 0, 0, 0);
	104	p1->SetBottomMargin(ratios ? 0.01 : .10);
	105	p1->SetFillColor(0);
	106	p1->SetBorderSize(0);
	107	p1->SetTopMargin(0.05);
	108	p1->SetRightMargin(0.05);
	109	p1->Draw();
	110	p1->cd();
0919aa56	111
	112	THStack* all = new THStack("all", "Analysis steps");
	113	if (res) {
8b1c306c	114	res->SetTitle("dN_{ch}/d#eta");
0919aa56	115	TIter next(res->GetHists());
	116	TH1* h = 0;
	117	while ((h = static_cast<TH1*>(next()))) all->Add(Rebin(h,rebin));
	118	}
	119	if (mcRes) {
8b1c306c	120	mcRes->SetTitle("Track-Refs");
0919aa56	121	TIter next(mcRes->GetHists());
	122	TH1* h = 0;
	123	while ((h = static_cast<TH1*>(next()))) all->Add(Rebin(h,rebin));
	124	}
	125	if (deltas) {
8b1c306c	126	deltas->SetTitle("#sum_{} #Delta/#Delta_{mip}");
0919aa56	127	TIter next(deltas->GetHists());
	128	TH1* h = 0;
	129	while ((h = static_cast<TH1*>(next()))) {
	130	h->SetMarkerStyle(25);
	131	all->Add(Rebin(h,rebin));
	132	}
	133	}
	134	if (nchs) {
8b1c306c	135	nchs->SetTitle("#sum_{} N_{ch,incl}");
0919aa56	136	TIter next(nchs->GetHists());
	137	TH1* h = 0;
	138	while ((h = static_cast<TH1*>(next()))) {
	139	h->SetMarkerStyle(21);
	140	all->Add(Rebin(h,rebin));
	141	}
	142	}
	143	if (prims) {
8b1c306c	144	prims->SetTitle("#sum_{} N_{ch,primary}");
0919aa56	145	TIter next(prims->GetHists());
	146	TH1* h = 0;
	147	while ((h = static_cast<TH1*>(next()))) {
	148	h->SetMarkerStyle(22);
	149	all->Add(Rebin(h,rebin));
	150	}
	151	}
	152	if (sumEta) all->Add(sumEta);
	153	all->Draw("nostack");
	154	all->GetHistogram()->SetXTitle("#eta");
	155	all->GetHistogram()->SetYTitle("signal");
	156	all->GetHistogram()->GetXaxis()->SetLabelFont(132);
	157	all->GetHistogram()->GetXaxis()->SetTitleFont(132);
	158	all->GetHistogram()->GetYaxis()->SetLabelFont(132);
	159	all->GetHistogram()->GetYaxis()->SetTitleFont(132);
	160	c->SetGridx();
	161
	162	TLegend* l = new TLegend(.33, .2, .53, .9);
8b1c306c	163	TLegendEntry* e = 0;
0919aa56	164	l->SetFillColor(0);
	165	l->SetFillStyle(0);
	166	l->SetBorderSize(0);
	167	l->SetNColumns(1);
	168	l->SetTextFont(132);
	169	Int_t i = 0;
	170	if (res) {
	171	TIter next(res->GetHists());
	172	TH1* h = 0;
	173	while ((h = static_cast<TH1*>(next()))) {
8b1c306c	174	e = l->AddEntry(Form("dummy%02d", i++),h->GetTitle(),"pl");
0919aa56	175	e->SetMarkerStyle(20);
	176	e->SetMarkerColor(h->GetMarkerColor());
	177	}
8b1c306c	178	e = l->AddEntry(Form("dummy%02d", i++),res->GetTitle(),"pl");
	179	e->SetMarkerStyle(20);
	180	e->SetMarkerColor(kBlack);
0919aa56	181	}
0919aa56	182	if (deltas) {
8b1c306c	183	e = l->AddEntry(Form("dummy%02d", i++), deltas->GetTitle(),"pl");
0919aa56	184	TH1* h = static_cast<TH1*>(deltas->GetHists()->At(0));
8b1c306c	185	e->SetMarkerStyle(h->GetMarkerStyle());
8b1c306c	186	e->SetMarkerColor(kBlack);
0919aa56	187	}
0919aa56	188	if (nchs) {
8b1c306c	189	e = l->AddEntry(Form("dummy%02d",i++),nchs->GetTitle(),"pl");
0919aa56	190	TH1* h = static_cast<TH1*>(nchs->GetHists()->At(0));
8b1c306c	191	e->SetMarkerStyle(h->GetMarkerStyle());
8b1c306c	192	e->SetMarkerColor(kBlack);
0919aa56	193	}
0919aa56	194	if (prims) {
8b1c306c	195	e = l->AddEntry(Form("dummy%02d", i++), prims->GetTitle(),"pl");
0919aa56	196	TH1* h = static_cast<TH1*>(prims->GetHists()->At(0));
8b1c306c	197	e->SetMarkerStyle(h->GetMarkerStyle());
8b1c306c	198	e->SetMarkerColor(kBlack);
0919aa56	199	}
	200
	201	if (mcRes) {
8b1c306c	202	e = l->AddEntry(Form("dummy%02d", i++), mcRes->GetTitle(), "pl");
0919aa56	203	TH1* h = static_cast<TH1*>(mcRes->GetHists()->At(0));
8b1c306c	204	e->SetMarkerStyle(h->GetMarkerStyle());
8b1c306c	205	e->SetMarkerColor(kBlack);
0919aa56	206	}
	207
	208	if (sumEta) l->AddEntry(sumEta);
	209	l->Draw();
	210
8b1c306c	211
	212	if (!ratios) return;
	213
	214	c->cd();
	215	TPad* p2 = new TPad("p2", "p2", 0, 0, 1, y1, 0, 0, 0);
	216	p2->SetTopMargin(0);
	217	p2->SetFillColor(0);
	218	p2->SetBorderSize(0);
	219	p2->SetRightMargin(0.05);
	220	p2->Draw();
	221	p2->cd();
	222
	223	THStack* rs = new THStack("ratios", "Ratios");
	224	Int_t nDN = Ratio(rs, deltas, nchs);
	225	Int_t nNR = Ratio(rs, nchs, res);
	226	Int_t nRP = Ratio(rs, res, prims);
	227	rs->Draw("nostack");
	228
	229	TLegend* ll = new TLegend(.38, .2, .48, .9);
	230	ll->SetFillColor(0);
	231	ll->SetFillStyle(0);
	232	ll->SetBorderSize(0);
	233	ll->SetNColumns(1);
	234	ll->SetTextFont(132);
	235	if (nDN) {
	236	e = ll->AddEntry("d1", Form("#frac{%s}{%s}",
	237	nchs->GetTitle(),
	238	deltas->GetTitle()), "lp");
	239	e->SetMarkerStyle(21);
	240	}
	241	if (nNR) {
	242	e = ll->AddEntry("d2", Form("#frac{%s}{%s}", res->GetTitle(),
	243	nchs->GetTitle()), "lp");
	244	e->SetMarkerStyle(20);
	245	}
	246	if (nRP) {
	247	e = ll->AddEntry("d3", Form("#frac{%s}{%s}", prims->GetTitle(),
	248	res->GetTitle()), "lp");
	249	e->SetMarkerStyle(22);
	250	}
	251	ll->Draw();
	252
	253
0919aa56	254	}
	255
	256
	257
	258