c287933545d8166505a1e263786eefc349150e53
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / scripts / FitELoss.C
1 #ifndef __CINT__
2 # include "AliForwardUtil.h"
3 # include <TFile.h>
4 # include <TList.h>
5 # include <TH1.h>
6 # include <TF1.h>
7 # include <TFitResult.h>
8 # include <TMath.h>
9 # include <TStyle.h>
10 # include <TArrow.h>
11 # include <TCanvas.h>
12 #else
13 class TCanvas;
14 class TFile;
15 class TH1;
16 class TList;
17 class TF1;
18 #endif
19
20 //____________________________________________________________________
21 TH1* GetEDist(TList* ef, UShort_t d, Char_t r, UShort_t etabin)
22 {
23   TList* dL = static_cast<TList*>(ef->FindObject(Form("FMD%d%c",d,r)));
24   if (!dL) {
25     Error("GetEDist", "Couldn't find list FMD%d%c",d,r);
26     ef->ls();
27     return 0;
28   }
29   if (etabin > 999) {
30     TH1* hist = static_cast<TH1*>(dL->FindObject(Form("FMD%d%c_edist",d,r)));
31     if (hist) {
32       Error("GetEDist", "Couldn't find EDists histogram for FMD%d%c",d,r);
33       return 0;
34     }
35   }
36       
37   TList* edL = static_cast<TList*>(dL->FindObject("EDists"));
38   if (!edL) {
39     Error("GetEDist", "Couldn't find list EDists for FMD%d%c",d,r);
40     return 0;
41   }
42   
43   TH1* hist = static_cast<TH1*>(edL->FindObject(Form("FMD%d%c_etabin%03d", 
44                                                      d, r, etabin)));
45   if (!hist) {
46     Error("GetEDist", "Couldn't find histogra FMD%d%c_etabin%03d",
47           d,r, etabin);
48     return 0;
49   }
50   
51   return hist;
52 }
53
54 //____________________________________________________________________
55 TH1* GetEDist(TList* ef, UShort_t d, Char_t r, Float_t eta)
56 {
57   if (!ef) { 
58     Error("GetEDist", "EF not set");
59     return 0;
60   }
61   TAxis* etaAxis = static_cast<TAxis*>(ef->FindObject("etaAxis"));
62   if (!etaAxis) { 
63     Error("GetEDist", "Couldn't find eta axis in list");
64     return 0;
65   }
66
67   UShort_t bin = etaAxis->FindBin(eta);
68   if (bin <= 0 || bin > etaAxis->GetNbins()) { 
69     Error("GetEDist", "eta=%f out of range [%f,%f] - getting ring histo", 
70           eta, etaAxis->GetXmin(), etaAxis->GetXmax());
71     return GetEDist(ef, d, r, UShort_t(1000));
72   }
73
74   return GetEDist(ef, d, r, bin);
75 }
76
77 //____________________________________________________________________
78 TList* GetEF(TFile* file) 
79 {
80   TList* forward = static_cast<TList*>(file->Get("PWG2forwardDnDeta/Forward"));
81   if (!forward) {
82     Error("GetEF", "Failed to get list PWG2forwardDnDeta/Forward from file");
83     return 0;
84   }
85   TList* ef = static_cast<TList*>(forward->FindObject("fmdEnergyFitter"));
86   if (!ef) {
87     Error("GetEF", "Failed to get energy fitter list");
88     return 0;
89   }
90   
91   return ef;
92 }
93
94 //____________________________________________________________________
95 TList* ef = 0;
96
97 //____________________________________________________________________
98 TList*  CheckEF()
99 {
100   if (ef) return ef;
101   
102   TFile* file = TFile::Open("AnalysisResults.root", "READ");
103   if (!file) { 
104     Error("Fit1", "Failed to open file");
105     return 0;
106   }
107   return GetEF(file);
108 }
109
110 //____________________________________________________________________
111 TCanvas* c = 0;
112
113 //____________________________________________________________________
114 TCanvas* CheckC()
115 {
116   if (c) return c;
117
118   gStyle->SetOptFit(1111);
119   gStyle->SetCanvasColor(0);
120   gStyle->SetCanvasBorderSize(0);
121   gStyle->SetCanvasBorderMode(0);
122   gStyle->SetCanvasDefW(800);
123   gStyle->SetCanvasDefH(800);
124   gStyle->SetPadTopMargin(0.05);
125   gStyle->SetPadRightMargin(0.05);
126   gStyle->SetTitleBorderSize(0);
127   gStyle->SetTitleFillColor(0);
128   gStyle->SetTitleStyle(0);
129   gStyle->SetStatBorderSize(1);
130   gStyle->SetStatColor(0);
131   gStyle->SetStatStyle(0);
132   gStyle->SetStatX(0.95);
133   gStyle->SetStatY(0.95);
134   gStyle->SetStatW(0.15);
135   gStyle->SetStatH(0.15);
136
137   c = new TCanvas("c", "c");
138   c->SetLogy();
139
140   return c;
141 }
142
143 //____________________________________________________________________
144 void PrintFit(TF1* f)
145 {
146   Int_t    nu     = f->GetNDF();
147   Double_t chi2   = f->GetChisquare();
148   Double_t chi2nu = (nu > 0 ? chi2/nu : 0);
149   Printf("%s: chi^2/nu=%f/%d=%f [%f,%f]", 
150          f->GetName(), chi2, nu, chi2nu,
151          f->GetXmin(), f->GetXmax());
152   for (Int_t i = 0; i < f->GetNpar(); i++) { 
153     Double_t v = f->GetParameter(i);
154     Double_t e = f->GetParError(i);
155     Double_t r = 100*(v == 0 ? 1 : e / v);
156     Printf("%32s = %14.7f +/- %14.7f (%5.1f%%)",f->GetParName(i),v,e,r);
157   }
158 }
159   
160 //____________________________________________________________________
161 void FitELoss(Int_t n, UShort_t d, Char_t r, Float_t eta)
162 {
163   TList* ef1 = CheckEF();
164   TCanvas* c1 = CheckC();
165   if (!ef1) return;
166   if (!c1)  return;
167
168   TH1* dist = GetEDist(ef1, d, r, eta);
169   if (!dist) return;
170
171   AliForwardUtil::ELossFitter f(0.4, 10, 4);
172   TF1* landau1 = new TF1(*f.Fit1Particle(dist, 0));
173   landau1->SetName("Landau1");
174   PrintFit(landau1);
175   TF1* landaun = new TF1(*f.FitNParticle(dist, n, 0)); 
176   landau1->SetName(Form("Landau%d", n));
177   PrintFit(landaun);
178   landau1->SetRange(0,10);
179   landaun->SetRange(0,10);
180   landau1->SetLineWidth(4);
181   landaun->SetLineWidth(4);
182   landau1->SetLineColor(kBlack);
183   landaun->SetLineColor(kBlack);
184
185   dist->GetListOfFunctions()->Clear();
186   dist->GetListOfFunctions()->Add(landau1);
187   dist->GetListOfFunctions()->Add(landaun);
188   dist->GetListOfFunctions()->ls();
189   dist->Draw();
190   landau1->Draw("same");
191   landaun->Draw("same");
192
193   Double_t mp = landaun->GetParameter(1);
194   Double_t xi = landaun->GetParameter(2);
195   for (Int_t i  = 1; i <= n; i++) { 
196     Double_t x  = i * (mp + xi * TMath::Log(i));
197     Double_t y  = landaun->Eval(x);
198     Double_t y1 = y < 0.05 ? 1 : 0.01;
199     TArrow* a = new TArrow(x,y1,x,y,0.03,"|>");
200     Info("FitSteer", "Delta_{p,%d}=%f", i, x);
201     a->SetLineWidth(2);
202     a->SetAngle(30);
203     a->Draw();
204   }
205
206   c1->cd();
207 }