d3276ac7dc6157018f8ef4d8998830ed791726f6
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / DrawFits.C
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 }