Added script to show fits
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / DrawFits.C
1 void
2 DrawFits(const char* fname="AnalysisResults.root")
3 {
4   TFile* file = TFile::Open(fname, "READ");
5   if (!file) {
6     Error("DrawFits", "Couldn't open %s", fname);
7     return;
8   }
9     
10   TList* forward = 
11     static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
12   if (!forward) { 
13     Error("DrawFits", "Couldn't get forward list from %s", fname);
14     return;
15   }
16
17   TList* fitter = 
18     static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
19   if (!fitter) { 
20     Error("DrawFits", "Couldn't get fitter folder");
21     return;
22   }
23   
24   TList stacks;
25   stacks.Add(fitter->FindObject("chi2"));
26   stacks.Add(fitter->FindObject("c"));
27   stacks.Add(fitter->FindObject("mpv"));
28   stacks.Add(fitter->FindObject("w"));
29   stacks.Add(fitter->FindObject("n"));
30   Int_t i=2;
31   while (true) { 
32     TObject* o = static_cast<THStack*>(fitter->FindObject(Form("a%d",i++)));
33     if (!o) break;
34     Info("DrawFits", "Adding %s", o->GetName());
35     stacks.Add(o);
36   }
37   // stacks.ls();
38   Int_t nMax = stacks.GetEntries();
39   for (Int_t i = nMax-1; i > 4; i--) { 
40     THStack* stack   = static_cast<THStack*>(stacks.At(i));
41     TIter    nextH(stack->GetHists());
42     TH1*     hist    = 0;
43     Bool_t   hasData = kFALSE;
44     while ((hist = static_cast<TH1*>(nextH()))) 
45       if (hist->Integral() > 0) hasData = kTRUE;
46     if (!hasData) nMax--;
47   }
48
49   gStyle->SetOptTitle(0);
50 #if 1
51   TCanvas* c = new TCanvas("c", "C",800,800);
52   c->SetFillColor(0);
53   c->SetBorderMode(0);
54   c->SetBorderSize(0);
55   c->Divide(1, nMax,0,0);
56
57   TIter next(&stacks);
58   THStack* stack = 0;
59   i = 1;
60   while ((stack = static_cast<THStack*>(next()))) {
61     if (i > nMax) break;
62     TVirtualPad* p = c->cd(i);
63     p->SetFillColor(0);
64     p->SetFillStyle(0);
65     p->SetGridx();
66     stack->Draw("nostack");
67     stack->GetHistogram()->SetYTitle(stack->GetTitle());
68     stack->GetHistogram()->SetXTitle("#eta");
69     TAxis* yaxis = stack->GetHistogram()->GetYaxis();
70     yaxis->SetTitleSize(0.15);
71     yaxis->SetLabelSize(0.08);
72     yaxis->SetTitleOffset(0.35);
73     yaxis->SetNdivisions(10);
74     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
75     xaxis->SetTitleSize(0.15);
76     xaxis->SetLabelSize(0.08);
77     xaxis->SetTitleOffset(0.35);
78     xaxis->SetNdivisions(320);
79     i++;
80     p->cd();
81   }
82 #endif
83     
84   gStyle->SetOptFit(111111);
85   gStyle->SetStatW(0.3);
86   gStyle->SetStatH(0.3);
87   gStyle->SetStatColor(0);
88   gStyle->SetStatBorderSize(1);
89   TCanvas* c1 = new TCanvas("c1", "c1", 800,800);
90   c1->SetFillColor(0);
91   c1->SetBorderMode(0);
92   c1->SetBorderSize(0);
93   c1->Divide(1, 5,0,0);
94
95   const char* dets[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
96   for (Int_t i = 0; i < 5; i++) { 
97     TVirtualPad* p = c1->cd(i+1);
98     p->SetGridx();
99     p->SetFillColor(0);
100     p->SetFillStyle(0);
101     TList* d = static_cast<TList*>(fitter->FindObject(dets[i]));
102     if (!d) { 
103       Warning("DrawFits", "List %s not found", dets[i]);
104       continue;
105     }
106     TH1* edist = static_cast<TH1*>(d->FindObject(Form("%s_edist", dets[i])));
107     if (!edist) {
108       Warning("DrawFits", "Histogram %s_edist not found", dets[i]);
109       continue;
110     }
111     edist->Draw();
112     TF1*   f = 0;
113     TIter  nextF(edist->GetListOfFunctions());
114     while ((f = static_cast<TF1*>(nextF()))) {
115       Double_t chi2 = f->GetChisquare();
116       Int_t    ndf  = f->GetNDF();
117       Printf("%s %s:\n  Range: %f-%f\n" 
118              "chi^2/ndf= %f / %d = %f", 
119              edist->GetName(), f->GetName(), 
120              f->GetXmin(), f->GetXmax(), chi2, ndf, 
121              (ndf > 0) ? chi2/ndf : 0);
122       for (Int_t j = 0; j < f->GetNpar(); j++) { 
123         Printf("  %-20s : %9.4f +/- %9.4f", 
124                f->GetParName(j), f->GetParameter(j), f->GetParError(j));
125       }
126     }
127     p->cd();
128   }
129   c1->cd();
130 }