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

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