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

#include <TFile.h>
#include <THStack.h>
#include <TList.h>
#include <TError.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TStyle.h>
#include <TF1.h>
#include <TLegend.h>
#include <TMath.h>

//____________________________________________________________________
// Global 
TList* fitter = 0;
TCanvas* canvas = 0;
const char* pdfName = "FitResults.pdf";

//____________________________________________________________________
TList* OpenFile(const char* fname)
{
  TFile* file = TFile::Open(fname, "READ");
  if (!file) {
    Error("DrawFits", "Couldn't open %s", fname);
    return 0;
  }
    
  TList* forward = 
    static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
  if (!forward) { 
    Error("DrawFits", "Couldn't get forward list from %s", fname);
    return 0;
  }

  fitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
  if (!fitter) { 
    Error("DrawFits", "Couldn't get fitter folder");
    return 0;
  }
  return fitter;
}
//____________________________________________________________________
TList* CheckFitter(const char* fname="AnalysisResults.root")
{
  if (!fitter) return OpenFile(fname);
  return fitter;
}

//____________________________________________________________________
TCanvas* CheckCanvas()
{
  if (canvas) return canvas;

  gStyle->SetOptFit(111111);
  gStyle->SetTitleFillColor(0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleX(0);
  gStyle->SetTitleY(.9);
  gStyle->SetTitleW(.4);
  gStyle->SetTitleH(.1);
  gStyle->SetStatW(0.2);
  gStyle->SetStatH(0.2);
  gStyle->SetStatColor(0);
  gStyle->SetStatBorderSize(1);
  gStyle->SetOptTitle(0);
  gStyle->SetOptFit(1111);
  gStyle->SetStatW(0.3);
  gStyle->SetStatH(0.15);
  gStyle->SetStatColor(0);
  gStyle->SetStatBorderSize(1);

  canvas = new TCanvas("c", "C", Int_t(800 / TMath::Sqrt(2)), 800);
  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);
  canvas->SetBorderSize(0);
  canvas->SetBottomMargin(0.15);
  return canvas;
}
  
//____________________________________________________________________
void DrawSummary(const char* fname="AnalysisResults.root")
{
  if (!CheckFitter(fname)) {
    Error("DrawFits", "File not opened");
    return;
  }
  if (!CheckCanvas()) {
    Error("DrawFits", "No canvas");
    return;
  }
  canvas->Clear();

  THStack* chi2nu;
  THStack* c;
  THStack* delta;
  THStack* xi;
  THStack* sigma;
  THStack* sigman;
  THStack* n;
  TList stacks;
  stacks.Add(chi2nu = static_cast<THStack*>(fitter->FindObject("chi2")));
  stacks.Add(c      = static_cast<THStack*>(fitter->FindObject("c")));
  stacks.Add(delta  = static_cast<THStack*>(fitter->FindObject("delta")));
  stacks.Add(xi     = static_cast<THStack*>(fitter->FindObject("xi")));
  stacks.Add(sigma  = static_cast<THStack*>(fitter->FindObject("sigma")));
  stacks.Add(sigman = static_cast<THStack*>(fitter->FindObject("sigman")));
  stacks.Add(n      = static_cast<THStack*>(fitter->FindObject("n")));
  Int_t baseA = stacks.GetEntries()+1;
  Int_t i=2;
  while (true) { 
    TObject* o = fitter->FindObject(Form("a%d",i++));
    if (!o) break;
    Info("DrawFits", "Adding %s", o->GetName());
    stacks.Add(o);
  }
  // stacks.ls();
  Int_t nMax = stacks.GetEntries();
  for (Int_t i = nMax-1; i >= baseA; i--) { 
    THStack* stack   = static_cast<THStack*>(stacks.At(i));
    TIter    nextH(stack->GetHists());
    TH1*     hist    = 0;
    Bool_t   hasData = kFALSE;
    while ((hist = static_cast<TH1*>(nextH()))) 
      if (hist->Integral() > 0) hasData = kTRUE;
    if (!hasData) nMax--;
  }

  canvas->SetRightMargin(0.05);
  canvas->SetTopMargin(0.05);
  canvas->Divide(2, (nMax+1)/2, 0.1, 0, 0);

  TIter next(&stacks);
  THStack* stack = 0;
  i = 0;
  Int_t b = 1;
  while ((stack = static_cast<THStack*>(next()))) {
    if (i > nMax) break;
    TVirtualPad* p = canvas->cd(1+i/5 + 2*(i%5));
    p->SetLeftMargin(.15);
    p->SetFillColor(0);
    p->SetFillStyle(0);
    p->SetGridx();
    stack->Draw("nostack");
    stack->GetHistogram()->SetYTitle(stack->GetTitle());
    stack->GetHistogram()->SetXTitle("#eta");

    TAxis* yaxis = stack->GetHistogram()->GetYaxis();
    if (i == 0)                     yaxis->SetRangeUser(0,20); // Chi2
    if (i == 1)                     stack->SetMaximum(1);      // c
    if (i == 2)                     stack->SetMaximum(1);      // delta
    if (i == 3)                     stack->SetMaximum(0.1);   // xi
    if (i == 4 || i == 5)           stack->SetMaximum(0.5);    // sigma{,n}
    if (i == 7)                     stack->SetMaximum(0.5);    // a
    yaxis->SetTitleSize(0.15);
    yaxis->SetLabelSize(0.08);
    yaxis->SetTitleOffset(0.35);
    yaxis->SetNdivisions(5);

    TAxis* xaxis = stack->GetHistogram()->GetXaxis();
    xaxis->SetTitleSize(0.15);
    xaxis->SetLabelSize(0.08);
    xaxis->SetTitleOffset(0.35);
    xaxis->SetNdivisions(10);

    // Redraw 
    stack->Draw("nostack");
    i++;
    if (i >= 5) b = 2;
    p->cd();
  }
  canvas->Print(pdfName, "Title:Fit summary");
}

//____________________________________________________________________
void DrawRings(const char* fname="AnalysisResults.root")
{
  if (!CheckFitter(fname)) {
    Error("DrawFits", "File not opened");
    return;
  }
  if (!CheckCanvas()) {
    Error("DrawFits", "No canvas");
    return;
  }
  canvas->Clear();
  
  canvas->Clear();
  canvas->Divide(1, 5,0,0);

  const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
  for (Int_t i = 0; i < 5; i++) { 
    TVirtualPad* p = canvas->cd(i+1);
    p->SetGridx();
    p->SetFillColor(0);
    p->SetFillStyle(0);
    TList* d = static_cast<TList*>(fitter->FindObject(dets[i]));
    if (!d) { 
      Warning("DrawFits", "List %s not found", dets[i]);
      continue;
    }
    TH1* edist = static_cast<TH1*>(d->FindObject(Form("%s_edist", dets[i])));
    if (!edist) {
      Warning("DrawFits", "Histogram %s_edist not found", dets[i]);
      continue;
    }
    edist->Draw();
    TF1*   f = 0;
    TIter  nextF(edist->GetListOfFunctions());
    while ((f = static_cast<TF1*>(nextF()))) {
      Double_t chi2 = f->GetChisquare();
      Int_t    ndf  = f->GetNDF();
      Printf("%s %s:\n  Range: %f-%f\n" 
             "chi^2/ndf= %f / %d = %f", 
	     edist->GetName(), f->GetName(), 
             f->GetXmin(), f->GetXmax(), chi2, ndf, 
	     (ndf > 0) ? chi2/ndf : 0);
      for (Int_t j = 0; j < f->GetNpar(); j++) { 
	Printf("  %-20s : %9.4f +/- %9.4f", 
	       f->GetParName(j), f->GetParameter(j), f->GetParError(j));
      }
    }
    p->cd();
  }
  canvas->cd();
  canvas->Print(pdfName, "Title:Fit to rings");
}

//____________________________________________________________________
void DrawEtaBins(const char* fname="AnalysisResults.root")
{
  if (!CheckFitter(fname)) {
    Error("DrawFits", "File not opened");
    return;
  }
  if (!CheckCanvas()) {
    Error("DrawFits", "No canvas");
    return;
  }
  canvas->Clear();
  canvas->Divide(2,2,0,0);

  for (UShort_t d=1; d<=3; d++) { 
    UShort_t nQ = (d == 1 ? 1 : 2);
    for (UShort_t q = 0; q < nQ; q++) { 
      Char_t r = (q == 0 ? 'I' : 'O');
      
      TList* ring = 
	static_cast<TList*>(fitter->FindObject(Form("FMD%d%c",d,r)));
      if (!ring) { 
	Error("PrintFits", "Couldn't get ring FMD%d%c", d,r);
	continue; 
      }
      TList* edists = static_cast<TList*>(ring->FindObject("EDists"));
      if (!edists) { 
	Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
	continue; 
      }
      TIter next(edists);
      TH1*  dist = 0;
      Int_t i    = 0;
      Int_t j    = 1;
      while ((dist = static_cast<TH1*>(next()))) { 
	if (i == 4) { 
	  i = 0;
	  j++;
	  canvas->Print(pdfName, Form("Title:FMD%d%c page %2d", d,r,j));
	}
	TVirtualPad* p = canvas->cd(++i);
	p->SetFillColor(kWhite);
	p->SetFillStyle(0);
	p->SetBorderSize(0);
	p->SetLogy();
	dist->SetMaximum(15);
	dist->Draw();
	
      }
      if (i != 0) {
	i++;
	for (; i <= 4; i++) {
	  TVirtualPad* p = canvas->cd(i);
	  p->Clear();
	  p->SetFillColor(kMagenta-3);
	  p->SetFillStyle(0);
	  p->SetBorderSize(0);
	}
	canvas->Print(pdfName, Form("FMD%d%c page %2d", d,r,j++));
      }
    }
  }
}   

//____________________________________________________________________
void
DrawFits(const char* fname="AnalysisResults.root")
{
  if (!CheckCanvas()) {
    Error("DrawFits", "No canvas");
    return;
  }
  canvas->Print(Form("%s[", pdfName));
  DrawSummary(fname);
  DrawRings(fname);
  DrawEtaBins(fname);
  canvas->Print(Form("%s]", pdfName));
}
Commit	Line	Data
c389303e	1	#include <TFile.h>
	2	#include <THStack.h>
	3	#include <TList.h>
	4	#include <TError.h>
	5	#include <TCanvas.h>
	6	#include <TPad.h>
	7	#include <TStyle.h>
	8	#include <TF1.h>
	9	#include <TLegend.h>
	10	#include <TMath.h>
	11
	12	//____________________________________________________________________
	13	// Global
	14	TList* fitter = 0;
	15	TCanvas* canvas = 0;
	16	const char* pdfName = "FitResults.pdf";
	17
	18	//____________________________________________________________________
	19	TList* OpenFile(const char* fname)
	20	{
	21	TFile* file = TFile::Open(fname, "READ");
	22	if (!file) {
	23	Error("DrawFits", "Couldn't open %s", fname);
	24	return 0;
	25	}
	26
	27	TList* forward =
	28	static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
	29	if (!forward) {
	30	Error("DrawFits", "Couldn't get forward list from %s", fname);
	31	return 0;
	32	}
	33
	34	fitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
	35	if (!fitter) {
	36	Error("DrawFits", "Couldn't get fitter folder");
	37	return 0;
	38	}
	39	return fitter;
	40	}
	41	//____________________________________________________________________
	42	TList* CheckFitter(const char* fname="AnalysisResults.root")
	43	{
	44	if (!fitter) return OpenFile(fname);
	45	return fitter;
	46	}
	47
	48	//____________________________________________________________________
	49	TCanvas* CheckCanvas()
	50	{
	51	if (canvas) return canvas;
	52
	53	gStyle->SetOptFit(111111);
	54	gStyle->SetTitleFillColor(0);
	55	gStyle->SetTitleBorderSize(0);
	56	gStyle->SetTitleX(0);
	57	gStyle->SetTitleY(.9);
	58	gStyle->SetTitleW(.4);
	59	gStyle->SetTitleH(.1);
	60	gStyle->SetStatW(0.2);
	61	gStyle->SetStatH(0.2);
	62	gStyle->SetStatColor(0);
	63	gStyle->SetStatBorderSize(1);
	64	gStyle->SetOptTitle(0);
65	gStyle->SetOptFit(1111);
66	gStyle->SetStatW(0.3);
67	gStyle->SetStatH(0.15);
68	gStyle->SetStatColor(0);
69	gStyle->SetStatBorderSize(1);
70
71	canvas = new TCanvas("c", "C", Int_t(800 / TMath::Sqrt(2)), 800);
72	canvas->SetFillColor(0);
73	canvas->SetBorderMode(0);
74	canvas->SetBorderSize(0);
75	canvas->SetBottomMargin(0.15);
76	return canvas;
77	}
78
79	//____________________________________________________________________
80	void DrawSummary(const char* fname="AnalysisResults.root")
81	{
82	if (!CheckFitter(fname)) {
83	Error("DrawFits", "File not opened");
84	return;
85	}
86	if (!CheckCanvas()) {
87	Error("DrawFits", "No canvas");
88	return;
89	}
90	canvas->Clear();
91
92	THStack* chi2nu;
93	THStack* c;
94	THStack* delta;
95	THStack* xi;
96	THStack* sigma;
97	THStack* sigman;
98	THStack* n;
99	TList stacks;
100	stacks.Add(chi2nu = static_cast<THStack*>(fitter->FindObject("chi2")));
101	stacks.Add(c = static_cast<THStack*>(fitter->FindObject("c")));
102	stacks.Add(delta = static_cast<THStack*>(fitter->FindObject("delta")));
103	stacks.Add(xi = static_cast<THStack*>(fitter->FindObject("xi")));
104	stacks.Add(sigma = static_cast<THStack*>(fitter->FindObject("sigma")));
105	stacks.Add(sigman = static_cast<THStack*>(fitter->FindObject("sigman")));
106	stacks.Add(n = static_cast<THStack*>(fitter->FindObject("n")));
107	Int_t baseA = stacks.GetEntries()+1;
108	Int_t i=2;
109	while (true) {
110	TObject* o = fitter->FindObject(Form("a%d",i++));
111	if (!o) break;
112	Info("DrawFits", "Adding %s", o->GetName());
113	stacks.Add(o);
114	}
115	// stacks.ls();
116	Int_t nMax = stacks.GetEntries();
117	for (Int_t i = nMax-1; i >= baseA; i--) {
118	THStack* stack = static_cast<THStack*>(stacks.At(i));
119	TIter nextH(stack->GetHists());
120	TH1* hist = 0;
121	Bool_t hasData = kFALSE;
122	while ((hist = static_cast<TH1*>(nextH())))
123	if (hist->Integral() > 0) hasData = kTRUE;
124	if (!hasData) nMax--;
125	}
126
127	canvas->SetRightMargin(0.05);
128	canvas->SetTopMargin(0.05);
129	canvas->Divide(2, (nMax+1)/2, 0.1, 0, 0);
130
131	TIter next(&stacks);
132	THStack* stack = 0;
133	i = 0;
134	Int_t b = 1;
135	while ((stack = static_cast<THStack*>(next()))) {
136	if (i > nMax) break;
137	TVirtualPad* p = canvas->cd(1+i/5 + 2*(i%5));
138	p->SetLeftMargin(.15);
139	p->SetFillColor(0);
140	p->SetFillStyle(0);
141	p->SetGridx();
142	stack->Draw("nostack");
143	stack->GetHistogram()->SetYTitle(stack->GetTitle());
144	stack->GetHistogram()->SetXTitle("#eta");
145
146	TAxis* yaxis = stack->GetHistogram()->GetYaxis();
147	if (i == 0) yaxis->SetRangeUser(0,20); // Chi2
148	if (i == 1) stack->SetMaximum(1); // c
149	if (i == 2) stack->SetMaximum(1); // delta
150	if (i == 3) stack->SetMaximum(0.1); // xi
151	if (i == 4 \|\| i == 5) stack->SetMaximum(0.5); // sigma{,n}
152	if (i == 7) stack->SetMaximum(0.5); // a
153	yaxis->SetTitleSize(0.15);
154	yaxis->SetLabelSize(0.08);
155	yaxis->SetTitleOffset(0.35);
156	yaxis->SetNdivisions(5);
157
158	TAxis* xaxis = stack->GetHistogram()->GetXaxis();
159	xaxis->SetTitleSize(0.15);
160	xaxis->SetLabelSize(0.08);
161	xaxis->SetTitleOffset(0.35);
162	xaxis->SetNdivisions(10);
163
164	// Redraw
165	stack->Draw("nostack");
166	i++;
167	if (i >= 5) b = 2;
168	p->cd();
169	}
170	canvas->Print(pdfName, "Title:Fit summary");
171	}
172
173	//____________________________________________________________________
174	void DrawRings(const char* fname="AnalysisResults.root")
175	{
176	if (!CheckFitter(fname)) {
177	Error("DrawFits", "File not opened");
178	return;
179	}
180	if (!CheckCanvas()) {
181	Error("DrawFits", "No canvas");
182	return;
183	}
184	canvas->Clear();
185
186	canvas->Clear();
187	canvas->Divide(1, 5,0,0);
188
189	const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
190	for (Int_t i = 0; i < 5; i++) {
191	TVirtualPad* p = canvas->cd(i+1);
192	p->SetGridx();
193	p->SetFillColor(0);
194	p->SetFillStyle(0);
195	TList* d = static_cast<TList*>(fitter->FindObject(dets[i]));
196	if (!d) {
197	Warning("DrawFits", "List %s not found", dets[i]);
198	continue;
199	}
200	TH1* edist = static_cast<TH1*>(d->FindObject(Form("%s_edist", dets[i])));
201	if (!edist) {
202	Warning("DrawFits", "Histogram %s_edist not found", dets[i]);
203	continue;
204	}
205	edist->Draw();
206	TF1* f = 0;
207	TIter nextF(edist->GetListOfFunctions());
208	while ((f = static_cast<TF1*>(nextF()))) {
209	Double_t chi2 = f->GetChisquare();
210	Int_t ndf = f->GetNDF();
211	Printf("%s %s:\n Range: %f-%f\n"
212	"chi^2/ndf= %f / %d = %f",
213	edist->GetName(), f->GetName(),
214	f->GetXmin(), f->GetXmax(), chi2, ndf,
215	(ndf > 0) ? chi2/ndf : 0);
216	for (Int_t j = 0; j < f->GetNpar(); j++) {
217	Printf(" %-20s : %9.4f +/- %9.4f",
218	f->GetParName(j), f->GetParameter(j), f->GetParError(j));
219	}
220	}
221	p->cd();
222	}
223	canvas->cd();
224	canvas->Print(pdfName, "Title:Fit to rings");
225	}
226
227	//____________________________________________________________________
228	void DrawEtaBins(const char* fname="AnalysisResults.root")
229	{
230	if (!CheckFitter(fname)) {
231	Error("DrawFits", "File not opened");
232	return;
233	}
234	if (!CheckCanvas()) {
235	Error("DrawFits", "No canvas");
236	return;
237	}
238	canvas->Clear();
239	canvas->Divide(2,2,0,0);
240
241	for (UShort_t d=1; d<=3; d++) {
242	UShort_t nQ = (d == 1 ? 1 : 2);
243	for (UShort_t q = 0; q < nQ; q++) {
244	Char_t r = (q == 0 ? 'I' : 'O');
245
246	TList* ring =
247	static_cast<TList*>(fitter->FindObject(Form("FMD%d%c",d,r)));
248	if (!ring) {
249	Error("PrintFits", "Couldn't get ring FMD%d%c", d,r);
250	continue;
251	}
252	TList* edists = static_cast<TList*>(ring->FindObject("EDists"));
253	if (!edists) {
254	Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
255	continue;
256	}
257	TIter next(edists);
258	TH1* dist = 0;
259	Int_t i = 0;
260	Int_t j = 1;
261	while ((dist = static_cast<TH1*>(next()))) {
262	if (i == 4) {
263	i = 0;
264	j++;
265	canvas->Print(pdfName, Form("Title:FMD%d%c page %2d", d,r,j));
266	}
267	TVirtualPad* p = canvas->cd(++i);
268	p->SetFillColor(kWhite);
269	p->SetFillStyle(0);
270	p->SetBorderSize(0);
271	p->SetLogy();
272	dist->SetMaximum(15);
273	dist->Draw();
274
275	}
276	if (i != 0) {
277	i++;
278	for (; i <= 4; i++) {
279	TVirtualPad* p = canvas->cd(i);
280	p->Clear();
281	p->SetFillColor(kMagenta-3);
282	p->SetFillStyle(0);
283	p->SetBorderSize(0);
284	}
285	canvas->Print(pdfName, Form("FMD%d%c page %2d", d,r,j++));
286	}
287	}
288	}
289	}
290
291	//____________________________________________________________________
292	void
293	DrawFits(const char* fname="AnalysisResults.root")
294	{
295	if (!CheckCanvas()) {
296	Error("DrawFits", "No canvas");
297	return;
298	}
299	canvas->Print(Form("%s[", pdfName));
300	DrawSummary(fname);
301	DrawRings(fname);
302	DrawEtaBins(fname);
303	canvas->Print(Form("%s]", pdfName));
304	}