55d242e3935759cd82c9a88a75f42a64d006f1a8
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / DrawFits.C
1 /**
2  * Script to draw the energy loss fits from the output file of 
3  * AliFMDELossFitter(Task). 
4  * 
5  * @ingroup pwg2_forward_analysis_scripts
6  */
7 #include <TFile.h>
8 #include <THStack.h>
9 #include <TList.h>
10 #include <TError.h>
11 #include <TCanvas.h>
12 #include <TPad.h>
13 #include <TStyle.h>
14 #include <TF1.h>
15 #include <TLegend.h>
16 #include <TMath.h>
17
18 //____________________________________________________________________
19 // Global 
20 TList* fitter = 0;
21 TCanvas* canvas = 0;
22 const char* pdfName = "FitResults.pdf";
23
24 //____________________________________________________________________
25 /** 
26  * Open a file.  The file is expected to contain the directory 
27  * structure 
28  *
29  * @verbatim 
30  *  file
31  *   +- Forward 
32  *       +- fmdEnergyFitter 
33  *           +- chi2   (THStack)
34  *           +- c      (THStack)
35  *           +- delta  (THStack)
36  *           +- xi     (THStack)
37  *           +- sigma  (THStack)
38  *           +- sigman (THStack)
39  *           +- n      (THStack)
40  *           +- a2     (THStack)
41  *           +- ...    (THStack)
42  *           +- an     (THStack)
43  *           +- FMD1I (TList)
44  *           |   +- FMD1I_edist (TH1)
45  *           |   +- EDists      (TList)
46  *           ...
47  * @endverbatim
48  * 
49  * @param fname File to open  
50  * 
51  * @return Pointer to the list of objects 
52  * 
53  * @ingroup pwg2_forward_analysis_scripts
54  */
55 TList* OpenFile(const char* fname)
56 {
57   TFile* file = TFile::Open(fname, "READ");
58   if (!file) {
59     Error("DrawFits", "Couldn't open %s", fname);
60     return 0;
61   }
62     
63   TList* forward = static_cast<TList*>(file->Get("Forward"));
64   // static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
65   if (!forward) { 
66     Error("DrawFits", "Couldn't get forward list from %s", fname);
67     return 0;
68   }
69
70   fitter = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
71   if (!fitter) { 
72     Error("DrawFits", "Couldn't get fitter folder");
73     return 0;
74   }
75   return fitter;
76 }
77 //____________________________________________________________________
78 /** 
79  * Open file if not done already 
80  * 
81  * @param fname File to open
82  * 
83  * @return List of fits 
84  * 
85  * @ingroup pwg2_forward_analysis_scripts
86  */
87 TList* CheckFitter(const char* fname="AnalysisResults.root")
88 {
89   if (!fitter) return OpenFile(fname);
90   return fitter;
91 }
92
93 //____________________________________________________________________
94 /** 
95  * Make canvas if not done already 
96  * 
97  * 
98  * @return Canvas 
99  * 
100  * @ingroup pwg2_forward_analysis_scripts
101  */
102 TCanvas* CheckCanvas()
103 {
104   if (canvas) return canvas;
105
106   gStyle->SetOptFit(111111);
107   gStyle->SetTitleFillColor(0);
108   gStyle->SetTitleBorderSize(0);
109   gStyle->SetTitleX(0);
110   gStyle->SetTitleY(.9);
111   gStyle->SetTitleW(.4);
112   gStyle->SetTitleH(.1);
113   gStyle->SetStatW(0.2);
114   gStyle->SetStatH(0.2);
115   gStyle->SetStatColor(0);
116   gStyle->SetStatBorderSize(1);
117   gStyle->SetOptTitle(0);
118   gStyle->SetOptFit(1111);
119   gStyle->SetStatW(0.3);
120   gStyle->SetStatH(0.15);
121   gStyle->SetStatColor(0);
122   gStyle->SetStatBorderSize(1);
123
124   canvas = new TCanvas("c", "C", Int_t(800 / TMath::Sqrt(2)), 800);
125   canvas->SetFillColor(0);
126   canvas->SetBorderMode(0);
127   canvas->SetBorderSize(0);
128   canvas->SetBottomMargin(0.15);
129   return canvas;
130 }
131   
132 //____________________________________________________________________
133 /** 
134  * Draw summary 
135  * 
136  * @param fname 
137  * 
138  * @ingroup pwg2_forward_analysis_scripts
139  */
140 void DrawSummary(const char* fname="AnalysisResults.root")
141 {
142   if (!CheckFitter(fname)) {
143     Error("DrawFits", "File not opened");
144     return;
145   }
146   if (!CheckCanvas()) {
147     Error("DrawFits", "No canvas");
148     return;
149   }
150   canvas->Clear();
151
152   THStack* chi2nu;
153   THStack* c;
154   THStack* delta;
155   THStack* xi;
156   THStack* sigma;
157   THStack* sigman;
158   THStack* n;
159   TList stacks;
160   stacks.Add(chi2nu = static_cast<THStack*>(fitter->FindObject("chi2")));
161   stacks.Add(c      = static_cast<THStack*>(fitter->FindObject("c")));
162   stacks.Add(delta  = static_cast<THStack*>(fitter->FindObject("delta")));
163   stacks.Add(xi     = static_cast<THStack*>(fitter->FindObject("xi")));
164   stacks.Add(sigma  = static_cast<THStack*>(fitter->FindObject("sigma")));
165   stacks.Add(sigman = static_cast<THStack*>(fitter->FindObject("sigman")));
166   stacks.Add(n      = static_cast<THStack*>(fitter->FindObject("n")));
167   Int_t baseA = stacks.GetEntries()+1;
168   Int_t i=2;
169   while (true) { 
170     TObject* o = fitter->FindObject(Form("a%d",i++));
171     if (!o) break;
172     Info("DrawFits", "Adding %s", o->GetName());
173     stacks.Add(o);
174   }
175   // stacks.ls();
176   Int_t nMax = stacks.GetEntries();
177   for (Int_t i = nMax-1; i >= baseA; i--) { 
178     THStack* stack   = static_cast<THStack*>(stacks.At(i));
179     TIter    nextH(stack->GetHists());
180     TH1*     hist    = 0;
181     Bool_t   hasData = kFALSE;
182     while ((hist = static_cast<TH1*>(nextH()))) 
183       if (hist->Integral() > 0) hasData = kTRUE;
184     if (!hasData) nMax--;
185   }
186
187   canvas->SetRightMargin(0.05);
188   canvas->SetTopMargin(0.05);
189   canvas->Divide(2, (nMax+1)/2, 0.1, 0, 0);
190
191   TIter next(&stacks);
192   THStack* stack = 0;
193   i = 0;
194   Int_t b = 1;
195   while ((stack = static_cast<THStack*>(next()))) {
196     if (i > nMax) break;
197     TVirtualPad* p = canvas->cd(1+i/5 + 2*(i%5));
198     p->SetLeftMargin(.15);
199     p->SetFillColor(0);
200     p->SetFillStyle(0);
201     p->SetGridx();
202     stack->Draw("nostack");
203     stack->GetHistogram()->SetYTitle(stack->GetTitle());
204     stack->GetHistogram()->SetXTitle("#eta");
205
206     TAxis* yaxis = stack->GetHistogram()->GetYaxis();
207     if (i == 0)                     yaxis->SetRangeUser(0,20); // Chi2
208     if (i == 1)                     stack->SetMaximum(1);      // c
209     if (i == 2)                     stack->SetMaximum(1);      // delta
210     if (i == 3)                     stack->SetMaximum(0.1);   // xi
211     if (i == 4 || i == 5)           stack->SetMaximum(0.5);    // sigma{,n}
212     if (i == 7)                     stack->SetMaximum(0.5);    // a
213     yaxis->SetTitleSize(0.15);
214     yaxis->SetLabelSize(0.08);
215     yaxis->SetTitleOffset(0.35);
216     yaxis->SetNdivisions(5);
217
218     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
219     xaxis->SetTitleSize(0.15);
220     xaxis->SetLabelSize(0.08);
221     xaxis->SetTitleOffset(0.35);
222     xaxis->SetNdivisions(10);
223
224     // Redraw 
225     stack->Draw("nostack");
226     i++;
227     if (i >= 5) b = 2;
228     p->cd();
229   }
230   canvas->Print(pdfName, "Title:Fit summary");
231 }
232
233 //____________________________________________________________________
234 /** 
235  * Draw ring fits 
236  * 
237  * @param fname 
238  * 
239  * @ingroup pwg2_forward_analysis_scripts
240  */
241 void DrawRings(const char* fname="AnalysisResults.root")
242 {
243   if (!CheckFitter(fname)) {
244     Error("DrawFits", "File not opened");
245     return;
246   }
247   if (!CheckCanvas()) {
248     Error("DrawFits", "No canvas");
249     return;
250   }
251   canvas->Clear();
252   
253   canvas->Clear();
254   canvas->Divide(1, 5,0,0);
255
256   const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
257   for (Int_t i = 0; i < 5; i++) { 
258     TVirtualPad* p = canvas->cd(i+1);
259     p->SetGridx();
260     p->SetFillColor(0);
261     p->SetFillStyle(0);
262     TList* d = static_cast<TList*>(fitter->FindObject(dets[i]));
263     if (!d) { 
264       Warning("DrawFits", "List %s not found", dets[i]);
265       continue;
266     }
267     TH1* edist = static_cast<TH1*>(d->FindObject(Form("%s_edist", dets[i])));
268     if (!edist) {
269       Warning("DrawFits", "Histogram %s_edist not found", dets[i]);
270       continue;
271     }
272     edist->Draw();
273     TF1*   f = 0;
274     TIter  nextF(edist->GetListOfFunctions());
275     while ((f = static_cast<TF1*>(nextF()))) {
276       Double_t chi2 = f->GetChisquare();
277       Int_t    ndf  = f->GetNDF();
278       Printf("%s %s:\n  Range: %f-%f\n" 
279              "chi^2/ndf= %f / %d = %f", 
280              edist->GetName(), f->GetName(), 
281              f->GetXmin(), f->GetXmax(), chi2, ndf, 
282              (ndf > 0) ? chi2/ndf : 0);
283       for (Int_t j = 0; j < f->GetNpar(); j++) { 
284         Printf("  %-20s : %9.4f +/- %9.4f", 
285                f->GetParName(j), f->GetParameter(j), f->GetParError(j));
286       }
287     }
288     p->cd();
289   }
290   canvas->cd();
291   canvas->Print(pdfName, "Title:Fit to rings");
292 }
293
294 //____________________________________________________________________
295 /** 
296  * Draw fits in eta bins 
297  * 
298  * @param fname 
299  * 
300  * @ingroup pwg2_forward_analysis_scripts
301  */
302 void DrawEtaBins(const char* fname="AnalysisResults.root")
303 {
304   if (!CheckFitter(fname)) {
305     Error("DrawFits", "File not opened");
306     return;
307   }
308   if (!CheckCanvas()) {
309     Error("DrawFits", "No canvas");
310     return;
311   }
312   canvas->Clear();
313   canvas->Divide(2,2,0,0);
314
315   for (UShort_t d=1; d<=3; d++) { 
316     UShort_t nQ = (d == 1 ? 1 : 2);
317     for (UShort_t q = 0; q < nQ; q++) { 
318       Char_t r = (q == 0 ? 'I' : 'O');
319       
320       TList* ring = 
321         static_cast<TList*>(fitter->FindObject(Form("FMD%d%c",d,r)));
322       if (!ring) { 
323         Error("PrintFits", "Couldn't get ring FMD%d%c", d,r);
324         continue; 
325       }
326       TList* edists = static_cast<TList*>(ring->FindObject("EDists"));
327       if (!edists) { 
328         Error("PrintFits", "Couldn't get EDists list for FMD%d%c", d,r);
329         continue; 
330       }
331       TIter next(edists);
332       TH1*  dist = 0;
333       Int_t i    = 0;
334       Int_t j    = 1;
335       while ((dist = static_cast<TH1*>(next()))) { 
336         if (i == 4) { 
337           i = 0;
338           j++;
339           canvas->Print(pdfName, Form("Title:FMD%d%c page %2d", d,r,j));
340         }
341         TVirtualPad* p = canvas->cd(++i);
342         p->SetFillColor(kWhite);
343         p->SetFillStyle(0);
344         p->SetBorderSize(0);
345         p->SetLogy();
346         dist->SetMaximum(15);
347         dist->Draw();
348         
349       }
350       if (i != 0) {
351         i++;
352         for (; i <= 4; i++) {
353           TVirtualPad* p = canvas->cd(i);
354           p->Clear();
355           p->SetFillColor(kMagenta-3);
356           p->SetFillStyle(0);
357           p->SetBorderSize(0);
358         }
359         canvas->Print(pdfName, Form("FMD%d%c page %2d", d,r,j++));
360       }
361     }
362   }
363 }   
364
365 //____________________________________________________________________
366 /** 
367  * Draw energy loss fits to a multi-page PDF
368  *
369  * The input file is the result of running AliFMDELossFitter - 
370  * either directly via AliFMDELossFitterTask or as part of a larger 
371  * train (AliForwardMultiplicityTask or AliForwardMCMultiplicityTask). 
372  * 
373  * @verbatim 
374  *  file
375  *   +- Forward 
376  *       +- fmdEnergyFitter 
377  *           +- chi2   (THStack)
378  *           +- c      (THStack)
379  *           +- delta  (THStack)
380  *           +- xi     (THStack)
381  *           +- sigma  (THStack)
382  *           +- sigman (THStack)
383  *           +- n      (THStack)
384  *           +- a2     (THStack)
385  *           +- ...    (THStack)
386  *           +- an     (THStack)
387  *           +- FMD1I (TList)
388  *           |   +- FMD1I_edist (TH1)
389  *           |   +- EDists      (TList)
390  *           ...
391  * @endverbatim
392  *
393  * @param fname 
394  * 
395  * @ingroup pwg2_forward_analysis_scripts
396  */
397 void
398 DrawFits(const char* fname="AnalysisResults.root")
399 {
400   if (!CheckCanvas()) {
401     Error("DrawFits", "No canvas");
402     return;
403   }
404   canvas->Print(Form("%s[", pdfName));
405   DrawSummary(fname);
406   DrawRings(fname);
407   DrawEtaBins(fname);
408   canvas->Print(Form("%s]", pdfName));
409 }
410 //
411 // EOF
412 //