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