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