Added some more scripts
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / DrawFits.C
1 /**
2  * Script to draw the energy loss fits 
3  * 
4  * @ingroup pwg2_forward_analysis_scripts
5  */
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     
32   TList* forward = static_cast<TList*>(file->Get("Forward"));
33   // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
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 }